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Abstract 



Einstein's system of equations in the ADM decomposition involves two subsystems of 
equations: evolution equations and constraint equations. For numerical relativity, one 
typically solves the constraint equations only on the initial time slice, and then uses the 
evolution equations to advance the solution in time. Our interest is in the case when 
the spatial domain is bounded and appropriate boundary conditions are imposed. A key 
difficulty, which we address in this thesis, is what boundary conditions to place at the 
artificial boundary that lead to long time stable numerical solutions. We develop an effective 
technique for finding well-posed constraint preserving boundary conditions for constrained 
first order symmetric hyperbolic systems. By using this technique, we study the preservation 
of constraints by some first order symmetric hyperbolic formulations of Einstein's equations 
derived from the ADM decomposition linearized around Minkowski spacctime with arbitrary 
lapse and shift perturbations, and the closely related question of their equivalence with the 
linearized ADM system. Our main result is the finding of well-posed maximal nonnegative 
constraint preserving boundary conditions for each of the first order symmetric hyperbolic 
formulations under investigation, for which the unique solution of the corresponding initial 
boundary value problem provides a solution to the linearized ADM system on polyhedral 
domains. 

We indicate how to transform first order symmetric hyperbolic systems with constraints 
into equivalent unconstrained first order symmetric hyperbolic systems (extended systems) 
by building-in the constraints. We analyze and prove the equivalence between the origi- 
nal and extended systems in both the case of pure Cauchy problem and initial boundary 
value problem. These results seem to be very useful for transforming constrained numeri- 
cal simulations into unconstrained ones. As applications, we derive the extended systems 
corresponding to the very same hyperbolic formulations of Einstein's equations for which 
boundary conditions consistent with the constraints have been found. Boundary conditions 
for these extended systems that make them equivalent to the original constrained systems 
are provided. 
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Chapter 1 



Introduction 

1.1 Background and Motivations 

In a nutshell, general relativity says that "matter tells space how to curve, while the cur- 
vature of space tells matter how to move," in a now famous phrase that the physicist John 
Wheeler once said. So, general relativity serves both as a theory of space and time and as a 
theory of gravitation. Einstein began with two basic but subtle and powerful ideas: gravity 
and acceleration are indistinguishable and matter in free-fall always takes the shortest pos- 
sible path in curved spacetime. One of Einstein's most important discoveries was a system 
of equations which relates spacetime and matter. 

While the theory of general relativity has tremendous philosophical implications and has 
given rise to exotic new physical concepts like black holes and dark matter, it is also crucial 
in some areas of modern technology such as global positioning systems. All these predic- 
tions and applications make general relativity a spectacularly successful theory. However, 
though a number of its major predictions have been carefully verified by experiments and 
observations, there are other key predictions, as the existence of gravitational waves, that 
remain to be fully tested. Einstein's equations possess solutions describing wavelike undula- 
tions in the spacetime. These gravitational waves correspond to ripples in spacetime itself, 
they are not waves of any substance or medium. Like electromagnetic waves, gravitational 
waves move at the speed of light and carry energy. In spite of carrying enormous amounts of 
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energy from some of the most violent events in the universe, gravitational waves are almost 
unobservable. For detecting and analyzing them, state-of-the-art detectors have been built 
in the United States and overseas. The development of these observatories, just coming on- 
line now, is one of the grandest scientific undertakings of our time, and the most expensive 
project ever funded by the National Science Foundation. With a network of gravitational 
waves detectors, mankind could open up a whole new window on the cosmic space. The 
study of the universe using gravitational waves would not be just a simple extension of the 
optical and electromagnetic possibilities, it would be the exploitation of an entirely new 
spectrum that could unveil parts and aspects of the universe inaccessible so far. 

This enormous technological effort to build ultra-sensitive detectors has been followed by 
an intense quest for developing computer methods to solve Einstein's equations. Having 
invested so much to detect gravitational waves, it is crucial that we be able to interpret the 
waveforms detected. Most recent investigations in numerical relativity have been based on 
first order hyperbolic formulations derived from the Arnowitt, Deser, and Misner, or ADM, 
decomposition 6, (see also j^) of Einstein's equations and some results have been obtained 
in spherical symmetry and axisymmetry. However, in the general three spatial dimensions 
case, which is needed for the simulation of realistic astrophysical systems, it has not been 
possible to obtain long term stable and accurate evolutions. One might argue that present 
day computational resources are still insufficient to carry out high enough resolution three 
dimensional simulations. However, the difficulty is likely to be more fundamental than that. 
It seems that there is insufficient understanding of the structure of Einstein's equations and 
there are too many unsolved questions related to how to approach them numerically. 

Einstein's system of equations can be decomposed into two subsystems of equations (ADM 
decomposition): evolution equations and constraint equations (Hamiltonian constraint and 
momentum constraints). For numerical relativity, one typically solves the constraint equa- 
tions only on the initial time slice, and then uses the evolution equations to advance the 
solution in time. A very difficult task is to derive good boundary conditions, and this 
problem is crucial if one takes into account that it seems impossible to have in the near 
future the computational power to put the boundaries far away from sources, far enough 
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that they would not affect the region of numerical spacetime being looked at. Traditionally, 
most numerical relativity treatments have been careful to impose initial data that satisfies 
the constraints. However, very rarely boundary conditions that lead to well-posedness are 
used and much less frequently are they consistent with the constraints. Stewart 47 has 
addressed this subject within Frittelli-Reula formulation linearized around flat space 
with unit lapse and zero shift in the quarter plane. Both main system and constraints 
propagate as first order strongly hyperbolic systems. This implies that vanishing values of 
the constraints at t = will propagate along characteristics. One wants the values of the 
incoming constraints at the boundary to vanish. However, one can not just impose them to 
vanish along the boundaries since the constraints involve derivatives of the fields across the 
boundary, not just the values of the fields themselves. If the Laplace-Fourier transforms are 
used, the linearity of the differential equations gives algebraic equations for the transforms 
of the fields. Stewart deduces boundary conditions for the main system in terms of Laplace- 
Fourier transforms that preserve the constraints by imposing the incoming modes for the 
system of constraints to vanish and translating these conditions in terms of Laplace-Fourier 
transforms of the main system variables. The idea of imposing the vanishing of the incom- 
ing constraints as boundary conditions is pursued further in ^1] within Einstein-Christoffel 
formulation [7j in the simple case of spherical symmetry. The radial derivative is eliminated 
in favor of time derivative in the expression of the incoming constraints by using the main 
evolution system. In these techniques are refined and employed for the linearized gen- 
eralized Einstein-Christoffel formulation [HH] around flat spacetime with vanishing lapse and 
shift perturbations on a cubic box. By considering well posed boundary conditions for the 
constraint system and trading normal derivatives for time and tangential ones, face systems 
are obtained that need to be solved first together with the compatibility conditions at the 
edges of the faces. The solutions of the face systems are used to impose well posed constraint 
preserving boundary conditions for the main system. A construction with several points 
in common with the one just described can be also found in 49 . These two papers, |15| 
and jlH]) are the closest to our work. A different approach can be found in j2,Sj . |24j . where 
the authors stray away from the general trend of seeking to impose the constraints along 
the boundary. Their method consists in making the four components of the Einstein tensor 
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projection along the normal to the boundary vanish. In the case of Einstein-Christoffel 
formulation restricted to spherical symmetry, the same boundary conditions as in are 
obtained. 

Before we end this brief review, it should also be mentioned here the work done on bound- 
ary conditions for Einstein's equations in harmonic coordinates |48j . |49j . when Einstein's 
equations become a system of second order hyperbolic equations for the metric components. 
The question of the constraints preservation does not appear here, as it is hidden in the 
gauge choice (the constraints have to be satisfied only at the initial surface, the harmonic 
gauge guarantees their preservation in time). 

What follows next is a summary of the contents of this dissertation, with emphasis on the 
ideas that connect the different parts. 

1.2 Thesis Organization 

This dissertation is divided into three main parts. The first part is represented by Chapter 
2 and mainly describes our results concerning first order symmetric hyperbolic (FOSH) sys- 
tems of partial differential equations. A special attention is being placed on FOSH systems 
with constraints and their well-posedness with or without boundary conditions. This first 
part represents a portion of the background theory needed for the rest of the dissertation. 
The second part. Chapter 3, is focused on the ADM decomposition of Einstein's equations 
due to Arnowitt, Deser, and Misner |2j and some important first order hyperbolic formula- 
tions derived from it. A novelty in this part is represented by the introduction and analysis 
in Subsection 3.4.4 of a new first order symmetric hyperbolic formulation of the linearized 
ADM decomposition due to Arnold The third part, represented by Chapter 4, is the 
most important part of this thesis. Here, we address a key difficulty in numerical relativity, 
the derivation of boundary conditions that lead to well posedness and consistent with the 
constraints. In the beginning of Chapter 4 we introduce and analyze a simpler model prob- 
lem which gives good insight for the more complex case of Einstein's equations. The core 
of Chapter 4 consists of the analysis of three important first order symmetric hyperbolic 
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formulations of Einstein's equations for which we provide well-posed constraint-preserving 
boundary conditions. 

In the remainder of this introduction, we will describe the principal results of this disserta- 
tion. 

1.3 Principal Results 

We have developed an effective and general technique for finding well-posed constraint 
preserving boundary conditions for constrained first order symmetric hyperbolic systems. 
The key point of this technique is the matching of the general forms of maximal nonnegative 
boundary conditions for the main system and the system of constraints. 

By applying this technique, we study the preservation of constraints by the linearized 
Einstein-Christoffel system around Minkowski spacetime with arbitrary lapse and shift per- 
turbations, and the closely related question of the equivalence of that system and the 
linearized ADM system. Our interest is in the case when the spatial domain is bounded 
and appropriate boundary conditions are imposed. However, we also consider the pure 
Cauchy problem with the result that the linearized Einstein-Christoffel and ADM systems 
are equivalent. Our main result is the finding of two distinct sets of well-posed maximal 
nonnegative constraint preserving boundary conditions for which the unique solution of 
the corresponding linearized Einstein-Christoffel initial boundary value problem provides a 
solution to the linearized ADM system on polyhedral domains. 

We have also obtained similar results for a very recent symmetric hyperbolic formulation 
of Einstein's equations introduced by Alekseenko and Arnold in 3 . A new first order 
symmetric hyperbolic formulation of linearized Einstein equations due to Arnold 2^ is an- 
alyzed. Again, the main result is the finding of well-posed constraint preserving boundary 
conditions. In fact, same ideas should be applicable to some other formulations and/or 
in different contexts, as, for example, linearization about some other backgrounds. How- 
ever, while the strategy of finding adequate boundary conditions is similar, the technical 
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apparatus employed depends very much on the formulation under investigation. 

Returning to the more general framework of constrained first order symmetric hyperbolic 
systems, we indicate how to transform such systems into equivalent unconstrained first 
order symmetric hyperbolic systems (extended systems) by building in the constraints. We 
also analyze and prove the equivalence between the original and extended systems in both 
the case of pure Cauchy problem and initial boundary value problem. As applications, we 
derive extended systems corresponding to the (EC) , (AA) , and (A) formulations respectively 
and boundary conditions that make them equivalent to the original constrained systems. 
These results seem to be useful for transforming constrained numerical simulations into 
unconstrained ones. 
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Chapter 2 



Symmetric Hyperbolic Systems 



2.1 Introduction 

In this chapter some basic results on first order symmetric hyperbohc (or FOSH) systems 
of partial differential equations are briefly reviewed, with special attention being given to 
systems with constraints and boundary conditions. All these results represent background 
material relevant to the discussions of the hyperbolic formulations of the Einstein equations 
which follow in the next chapters. Much more information on hyperbolic systems can be 
found in the books by John [SSI, Kreiss and Lorenz [SS], Gustafsson, Kreiss and Oliger fI7\ . 
and Evans jJTj, among many others. 

The second section of this chapter is intended to enlist the basic definition of FOSH systems 
and some relevant existence and uniqueness results. The third section is dedicated to the 
analysis of constrained initial value problems in a more abstract framework and in the case of 
FOSH systems of partial differential equations. The emphasis is on the equivalence between 
a given system subject to constraints and a corresponding extended unconstrained system. 
The fourth section deals with boundary conditions for FOSH systems and the connections 
between the initial boundary value problem for a given FOSH system with constraints and 
that for the extended system. Section 2.3 and a substantial part of Section 2.4 represent 
our contribution to the subject. 
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2.2 Initial Value Problems 



In this section we will be concerned with a linear first order system of equations for a column 
vector u = u{x, t) = u{xi, . . . ,X]\f,t) with m components ui, . . . , Um- Such a system can be 
written as 

TV 

Lu = dtu + ^Ai{x,t)diU = f{x,t) in x (0,T], (2.1) 

i=l 

where T > 0. Here Ai, . . . , A]\r are given m x m matrix functions, and / is a given m- 
vector field. We will further assume that A^ are of class C^, with bounded derivatives over 
X [0,r], and / G H\R^ x (0,r);M™). 

As initial data we prescribe the values of u on the hyperplane t = 

It = no on X {t = 0}, (2.2) 
with no e i7^(M^;M™). For each n; e M^, define 

N 

A{w){x, t) =Y^ WiAi{x, t) {x e M^, t > 0). (2.3) 

i=l 

The system ()2.1() is called symmetric hyperbolic if Ai[x^t) is a symmetric m x m matrix 
for each x G M^, t > [i = 1,...,N). Thus, the m x m matrix A[w){x^t) has only 
real eigenvalues and the corresponding eigenvectors form a basis of for each w G M^, 
X G M^, and t > 0. 

Remark 1. More general systems having the form 

N 

Aodtu + Ai{x, t)diu + B{x, t)u = f{x, t) (2.4) 

i=l 

are also called symmetric, provided the matrix functions Ai are symmetric for i = 0, . . . , N , 
and Aq is positive definite. The results set forth below can be easily extended to such systems. 
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As in ilTj, Section 7.3.2., let us define the bilinear form 



A[u,v;t] := / > {Ai{- , t)diu) ■ v dx 



for < t < T, n, ?; G H^{M.^-W^). 



Definition. We say 



u G L2((0,T); i/i(M^;E™)), with u' G L2((0,T); L2(M^;M™)) 



is a weak solution of the initial value problem (|2.1|) . (|2.2j) provided 

(i) (n',^) +^['u,t;;t] = {f,v) for each G H'^{M.^;W) and a.e. < t < T, and 

(ii) 'u(O) = uq. 

Here ( • , • ) denotes the inner product in L2(M^;M™). 

By using energy methods and the vanishing viscosity technique (see |17|. Section 7.3.2.), 
the following existence and uniqueness result can be proven: 

Theorem 1. The initial value problem (|2.1j) . (j2.2|) has a unique weak solution. 

In what follows, we will be more interested in first order symmetric hyperbolic systems with 
constant coefficients. For such systems, a more general result (including regularity) is valid. 

Theorem 2. /i?/, Section 1.3.3. (also JM, Section 16.1.) AssumeuQ G i7"(M^;M'"), with 
s > N/2 + m. Then there is a unique solution u G C"'^(]R^ x [0, oo); M™') of the initial value 
problem (jHU), (1^ . 

The main tool used for proving this theorem is the Fourier transform. The unique 
solution is given by: 



ti(x, t) 



(27r)^/2 J^M 



[ e"-'"[e^*^('")uo(?^) + [ e- 



;^(*-^)^(^)/(^l;,s)ds] dw. 



(2.5) 
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2.3 Constrained Initial Value Problems 



2.3.1 Abstract Framework 

We introduce in this subsection an extended system corresponding to a given constrained 
system defined on Hilbert spaces and investigate the equivalence of these two systems. 

Let us consider the following system subject to constraints: 

y = Ay + f, (2.6) 
By = 0, (2.7) 
2/(0) = 2/0, (2.8) 

where A : D[A) C y ^ y, B : D{B) CY^X are densely defined closed linear operators 
on the Hilbert spaces {X, ( • , • )x) and {Y, {■ , ■ )y), and / : [0, oo) Y. Moreover, suppose 
A is skew-symmetric and 

y4(Ker5) C Ker5. (2.9) 

Of course, we assume that the compatibility condition ByQ = is satisfied. Moreover, 
another more subtle compatibility condition must hold: Bf = 0. This is because, for any 
fixed t G [0, oo), B{[y{t) — y{t)]/{t — t)) = 0, for all t > 0, and passing to the limit as t — > 
it turns out that B{y{t)) =0, Vt € [0, oo) (since i? is a closed operator). By operating on 
((T^ with B, it follows that 

= B{y{t)) = BAy{t) + Bf{I) = Bf{t), Vt G [0, oo), 

where the last equality comes from ()2.9|) . 

Remark 2. If f = 0, then the energy of the solution is preserved: 

m = l\\yit)f = l\\yof, vt>o. 
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Theorem 3. The pair (y, z)^ is a solution of the associated unconstrained system 

:)=(;■:)(■)•(:) 

y(0) = yo, z{0) = 0, (2.11) 
if and only if z = and y satisfies the constrained system (|2.6j) - (|2.8|) . 

Proof If z = and y satisfies (ESJ-jlSl), then it is clear that (y,0) satisfies (piIH) - (|TTT]) . 

Now suppose {y,z) satisfies (|2.1fl|) ~ (j2.1Hl . Observe that we can spUt y as a direct sum as 
fohowing 



Y = KeTB e (KerB)^ = KerB eimB*. (2.12) 

According to this decomposition of Y, 

y{t) = yi{t) + y2{t), (2.13) 



where yi{t) G Keri?, and 2/2 (*) £ liaB*. 



Since both Ker B and Im B* are closed and the corresponding projections are continuous, 

y(t) = yi(t) + y2(t), (2.14) 



with yi{t) G KerS, and 2/2 (t) G ImB*. 

Prom (HHni), and (fTTT|) . we obtain that 



yi = Ayi + f, (2.15) 
yi(0)=yo, (2.16) 
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and 



m = Ay2 - B*z, 



(2.17) 



z = By2 



(2.18) 



2/2(0) = 0, z(0) = 0. 



(2.19) 



Observe that 





Therefore, 



z{t)f + \\y2{t)r 



^(0)f + ||2/2(0)f =0, 



which imphes z = and 7/2 = 0. Thus, y = yi and (|2.6|) - ()2.8() are satisfied. 



□ 



2.3.2 Constrained First Order Symmetric Hyperbolic Problems 

In this subsection, we will prove a result similar to Theorem |31 for the initial value problem 



where A = "^^^i Aidi, with Ai £ ]^'T«xm constant symmetric matrices, and B = "^^^i Bidi, 
with Bi E RP^"* constant matrices. Of course, we assume that 1)2. 9() and the compatibility 
conditions Buq = 0, Bf{-, t) = 0, Vt > 0, hold. 

Theorem 4. (equivalence for classical solutions) IJuq G F'*(IR^; W^) andf{-,t) G F'*(M^; M™), 
yt > 0, for s > N/2 + m, then the pair {u,z f G C^{R^ x [0, 00); M'"+p) is a solution of 



dtu = Au + f, 



(2.20) 



Bu = 



(2.21) 



u{x, 0) = no 



(2.22) 
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the associated unconstrained system 



d u \ A -B* \ u \ / / , 

' * ' * ' * + r , (2-23) 



dt 



B J \ z I \ 

u{x,0) = uo, z{x,0) =0, (2.24) 



if and only if z = 0, and u E C"'^(M^ x [0,00);!^"^) satisfies the constrained system ()2.2U() - 

Proof If z = and -u satisfies (|2?2n|) - (fz!22|) . then it is clear ttiat (y,0)'^ satisfies (j22S|- 
Now, let us prove the converse. Denote by 

\b J fri dx, 

From ()2.5p . we know that the solution of (|2.23|) - H2.24() has the following expression 



^ = /m^v e-ne-*^(^) I ]{y) 

_^jt^-i(ts)My) I I {y,s)ds]dy. 



(2.25) 



The next step in the proof is to show that 



Ai,) r; = ^"^'j'^ I (2.26) 



for all positive integer n. 



We are going to prove (|2.26|) by induction 
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For n = 1, we have 

A(y) f "° 1 = f 

But, since Buq = 0, it follows that B(y)uQ = by taking the Fourier transform. So 



My) 



^0 \ _ I A{y)'^o 
j I 



Assume that (|2.26|) is true for n = k — 1 and let us prove it for n = k. 

\ J \ J \ B(!,)/l'--i(!;)«o 

Since uq e Ker5, from (|2.9|) . we can see that 

BA''~\o = 0. (2.27) 

Applying the Fourier transform to (|2.27j) . we get 

Biy)A''-\y)uo = 0. 

Thus, 



and the proof of 1)2. 26(1 is complete. 



/ \ 



From (|2.26|) . observe that 



n! 



n=0 



/ no 1 






I ) 
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Same arguments show that 



-i{t-s)A{y) 



/ \ / V 



From (lOSll and 12123), it follows that 



(; 


i-i 











with 

"(^' = lArm I e-^[e-*-^(j')no(y) + /* e-(*-^)^(^)/(y, s)ds] dy 
[I-kY"/^ Jm.n Jo 

the unique solution of (tT^ - (IT^ . □ 



2.4 Boundary Conditions 

In general, one has to be careful when choosing boundary conditions for a hyperbolic equa- 
tion (or system). This can be seen even in the simple case of a first order equation in one 
space dimension (the transport equation). It seems that any acceptable boundary condi- 
tions should give the incoming modes into the spatial domain, but they must not try to 
change the behavior of the outgoing modes. In several dimensions the situation is much 
more complicated since there is not easy to identify the incoming and outgoing modes. 
Worse, there may also be waves moving tangent to the boundary, and it is not very clear 
how these modes could be casted into the boundary conditions. 

A few approaches to these questions have been proposed. Some answers have been given by 
Friedrichs ^H] via the "energy method" (see also the work done by Courant and Hilbert). 
This method provides criteria which are sufficient for constructing boundary conditions 
that lead to a well-posed problem. Other sufficient conditions have been pointed out by 
Lax and Phillips in their very interesting work jlS7| . A necessary and sufficient condition 
for having a well-posed initial boundary condition has been proved by Hersh (22j, but his 
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result was only for systems with constant coefficients and defined on a half-space with non- 
characteristic boundary conditions. Using Fourier and Laplace transforms, he constructed 
solutions and derived a necessary and sufficient condition for well-posedness. Later on 
(in the 1970s and 1980s), more technical approaches came up. Kreiss [34 , Majda and 
Osher [32], among others, proved pretty complicated algebraic results concerning boundary 
conditions. Remarkably, Kreiss [25 gave a criteria that determine whether a boundary 
condition is admissible or not. The main point of his approach was the possibility to 
solve for incoming modes in terms of outgoing modes and boundary conditions. Majda 
and Osher [HHl generalized Kreiss' theory to the case of uniformly characteristic boundary. 
Other significant contributions to this subject have been made by Rauch [301) Higdon |30j . 
Secchi |42j-|46j. among many others. 

2.4.1 Mciximal Non-Negative Boundary Conditions 

In this subsection we prove a well-posedness result for first order symmetric hyperbolic 
initial boundary value problems that closely follows the ideas of |37) . |40j . and jJH]- More- 
over, we give an algebraic characterization of maximal non-negative boundary conditions 
which will be used later for determining constraint preserving boundary conditions for some 
constrained first order symmetric hyperbolic systems. 

Consider the symmetric hyperbolic system of equations (|2.H) on 0,t = ^ x (0,T), where 
C is a bounded domain with a smooth boundary dil, and / G L^(r22^, R*"). 

Set n{x) = {rii, . . . ,nj\f) be the outer normal to at x G dO,, and denote by An{x,t) the 
boundary matrix 

N 

An{x,t) = Y,n^A\ (2.30) 
1=1 

We supplement 1)2.1(1 with the initial condition (|2.2|) . with uq G H^{Q), and with linear 
boundary conditions of the following form 

E{x, t)u{x, t) = 0, on dfl x (0, T). (2.31) 
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Of course, we suppose that the compatibihty condition E{x,0)uo{x) = holds on dQ. 

In fact, by choosing a function of //^(fij-, M™") satisfying (|2.2j) and (|2.31|) . and changing the 
variable, we may assume that uq = 0; hereafter, we stick with this choice of uq. 

Also, the boundary condition ()2.31|) may be regarded as 

u{x,t) E N{x,t) = Ker E{x,t), V(x,t) e dQ x {0,T). (2.32) 

Denote the formal adjoint of L by L* 

N 

L*u = -dtu-'^diiA^u). 

i=l 

Associated to ()2.1() . (|2.2() . and H2.31() (or H2.32() ). we consider the adjoint problem 

L*v = f innr, (2.33) 
v{x,t) e {An{x,t)N{x,t))-^, y{x,t) £dnx (0,T), (2.34) 
v{x,T) = inn. (2.35) 

Next, define the admissible spaces of solutions for both the original problem and the adjoint 
problem: 

n = {ue H\nT,'^"') ■■ u satisfies (j22I) and (tOTT) (or ^M))} 

and 

n^ = {v e H\nT,M'^) : V satisfies (fO^ and (tO^l) }. 

Observe that, if u G W and v £ Tl^:, then from Green's formula 

{v, Lu) — {u, L*v) = / Anuda, 
Jdnx(0,T) 
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it follows that 



{v, Lu) = (n, L*v). 



Definition. The function u G L^(Oj',M'") is said to be a weak solution of (|2.H) . (|2.2|) . and 
(tr^TTl) if 

(z;,/)-(L*7;,'u) = 0, Vf G W*. 

Definition. The function u £ L^(r2T,IR™') is said to be a strong solution of (|2.H) . (|2.2j) . 
and ()2.31|) . if the pair (n, /) belongs to the closure of the graph of L; in other words, if u is 
the limit in the norm of a sequence of functions {u^} C Ti such that = Lu^ ^ / in 
the L'^ norm. 

Remark 3. If u solves (|2.H) . (|2.2j> . and (|2.Hljl in i/ie strong sense, then it also solves the 
problem in the weak sense. 

Definition. The boundary condition (|2.31() (or H2.32() ) is called non-negative if the matrix 
An{x,t) is non-negative over N{x,t) 

-J AnU > 0, Vn G A^(x,t), G x (0,r). (2.36) 

Theorem 5. 14(J\I ) If the boundary condition ()2.31|) is maximal non-negative at each 

point {x,t) G x (0,T) (meaning that (|2.36|) holds and there is no other larger subspace 
containing N{x,t) and having the same property), then ^2.1^ . (j2.2j) . and (|2.3H) has a unique 
strong solution for any given integrable function f . 

In order to prove this theorem, we need a couple of intermediate results. 

Lemma 6. For all u £ H satisfying the condition (|2.36j) . there exists a positive constant 
C > that does not depend on u such that 

\\u\\L'2{nT,K"') ^ C\\Lu\\L2(^Q^^^n^y (2.37) 
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Proof. Symmetric operators with smooth satisfy the following identity 



I d , rr ^ 1 A a 



^^^^ = ^i^^"^^) + ^ ^ ^^^^^^^^ + (^-^^^ 

2=1 * 

where 

2=1 

If we integrate 1)2.38^ over = x (0, s), we get 

/ u^Lu = 11 u^u +]- f u^AnU + [ u^Ku. (2.40) 

JvLs ^ Jnx{s} ^ Jdnx{o,s) J^s 

From ()2.36p and H2.4U() . it is easy to see that 

[ u^u<2 [ u'^Lu + 2\\K\\ [ u^u<{l + 2\\K\\) [ u^u+ [ {Luf{Lu). (2.41) 

Denote by 

(j){s) = / u^u. 

Then ()2.4ip recasts into 

cl>\s) < (1 + 2\\Kms) + \\Lu\\l,^^^^^^y (2.42) 
Applying the Gronwall's lemma to ()2.42p . it follows that 

gS(l+2||E-||) 

<^(«) ^ i + 2||i^|| ll-^^ll^^(^T,M-)- 



Hence, 

^ ^T(l+2\\K\\) ^ 

(t){T) = ||'u|l2,2(nj,,K™) < 1 + 2||Er|| "^""-f^^(f^T,iR™)- 
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Then, the inequality H2.37() holds for 



C 



\^T{\+2\\K\\ 
l + 2||i^| 



□ 



A simple consequence of this lemma is stated next. 

Corollary 7. For any square integrahle function f there is at most one solution u & TC 
satisfying non-negative boundary conditions to the problem H2.1() . ()2.2|) . and H2.31() . 

Lemma 8. // the boundary conditions ()2.32|) are maximal non-negative for L, then the 
adjoint boundary conditions (|2.34() are non-negative for L* . 

Proof. According to the expression of the boundary matrix for is equal to — ^ri* 
Hence it suffices to show that v'^AnV is non-positive for each v satisfying the adjoint bound- 
ary condition. 

Arguing by contradiction, suppose that there is a v such that v'^AnV > 0. 

Consider the linear space N{x,t) © v. The elements in this space have the form u-\- av, 
where u G N{x, t) and a is a real number. 

Observe that 

(n -|- av)'^An{u + av) = u^AnU + 2av'^AnU + a'^v'^AnV = vF A^u -\- a'^v'^AnV > 0, 
which is in contradiction with the maximality of N{x,t). □ 

Proof, (of Theorem [SJ Existence. We claim that L{Ti.) is dense in L^(r27^,M™). Arguing 
by contradiction, let us suppose the contrary. Then there exists a non-trivial function v 
orthogonal to L(Ti.), i.e. 

{Lu,v) = 0, Vu G n. 
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Therefore, u is a weak solution to the adjoint problem corresponding to / = 0. 

Now, let us prove that v = 0. From Theorem 1.1, we get that f is a strong solution; 
therefore there exists a sequence f„, — > u in L^, such that /„ = Lvn — > in and Vn 
satisfies the (non-negative) adjoint boundary conditions, which imphes 

ll^'nll <C||M| ^0. 

So, v = 0. In conclusion, if / is any element of L^^Qt,^"^), then we can find a sequence 
{un} C H so that Lun ^ /, as n ^ oo, in L^. 

From Lemma IHl it follows that {n„} is a Cauchy sequence in L^. Thus, there is u E 
L^(i77-,M'" such that ^ u and Lun — > / in L^. 

This imphes that u is a strong solution of ((SIJ, (|T2|) . and (|OT|) for the given / e L?{Q.t, K™). 
Uniqueness and Continuous Dependence. Follow from Lemma IHl □ 

There are results concerning the regularity of the solution of ((211, (EIS), and (lOTjl . If 

/ G Li((0,r);L2(17)) and uq € L'^{^) then the solution u belongs to C((0, T); ^2(17)) (see 
|4()j). In fact, this regularity result can be improved by imposing more regularity for / and 
uq and, in addition, compatibility conditions at the corner {t = {)}xdVt. These compatibility 
conditions are computed in the following fashion. Denote by ■n{x) the orthogonal projection 
onto N{x)-^. The compatibility condition of order j comes from expressing dl{iTu) at 
{t = 0} X dQ in terms of uq and / and requiring that the resulting expression vanishes. For 
example, for j = 0, we find ttuq = on dQ, or uq G N{x) for all x S dQ. For j = 1, we 
have to impose: vtuq = and vr(/(0, •) — YliLi ^i^i'^o) = on dO,. 

Definition. ([10]) A smooth vector field 7 on is called tangential if and only if, for 
every x £ (7(x),n(x)) = 0. For s a positive integer, the space H^^^{Q) consists 
of those u E L2(J7) with the property that for any I < s and tangential fields {71}'=^, 
7i72 • • • 7zii G L'^{0,) (see for more on -ff^n spaces). 

Theorem 9. (I40\l ) Suppose s > 1 is an integer and Ai, N, 90 are of class C^'^. Suppose 
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the data uq £ and &[ f E L^{{0,T); H^^^ (Q)) for < j < s and in addition there is 
a < T' < T such that d{f G L^{{0,T'); H^~^n)), < j < s. If the data satisfy the 
compatibility conditions up to order s — 1, then the solution u of (|2.H) . (|2.2j) . and (|2.31|) 
Ues m C-'GO, T\-Hl~Jm for 0<j<s. 

We close this subsection by giving an algebraic characterization of maximal non-negative 
boundary conditions. 

Suppose that the boundary matrix An{x,t) has Iq 0-eigenvalues Xi{x,t),... ,Xii^{x,t), l- 
negative eigenvalues Xig+i{x,t), . . . ,Xig^i_{x,t), and /+ positive eigenvalues AjQ_|_i__|_i(x, t), 
. . . , Xm{x,t). Let ei{x,t), . . . , ei^, c/q+i, . . . , eig+i_, ej^+u+i, . . . , emix,t) be the corre- 
sponding eigenvectors. Naturally, at {x, t) £ dQ x (0, T), any vector of N{x, t) has the form 
V = YlT^i Oii^i{x^t)- The non-negative condition 1)2. 36() implies: 

m 

An{x, t)v = XiO^ > 0, 
i=l 

or 

^i^l > - E ^^"^ (2.43) 

Ai>0 Aj<0 

Observe that the dimension of A^(x, t) must be equal to the number of positive and null 
eigenvalues counted with their multiplicities. 

Now, any v G N{x,t) can be written as v = Qa, where a = {ai, . . . ,am)'^ , and Q = 
(ei(x,t), . . . ,em{x,t)). 

Since N{x, t) is a subspace of codimension there exists al^ xm matrix E{x, t) such that 

N{x,t) = {v : E{x,t)v = 0}. 

So, EQa = 0, or 

SoOq + S^a- — S'+a+ = 0, (2.44) 

where ao = (ai, . . . ,aZ(,), «_ = {ai^+i, . . . ,aig+i_), a+ = {aig+i_+i, . . . ,am), and Sq = 
EQo, S- = EQ^, S+ = -EQ+, with Qq = (ei(x, t), . . . ,eio(x, t)), Q_ = {ei^+i{x,t), . . . , 
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eig+i_{x,t)), and Q+ = {eiQ+i_+i{x,t), . . . ,em{x,t)). 

Since N{x,t) is maximal non-negative, it follows that KerAn{x,t) C N{x,t), and so, 



So = EQq = Oi_xio- 



Therefore, (|2.44|) reads 



— S'+a-L = 0. 



(2.45) 



Observe that the /_ x /_ matrix S- is invertible. If not, there exists a_ 7^ so that 
S-O- = 0, and so. 



lo+l- 

t> = ^ aiei{x,t) 

i=lo + l 



belongs to N(x,t). But this is in contradiction with (|2.43j) . Hence, S- is invertible. 



So, l|2.45|l recasts into 



Now, we must have ()2.43() satisfied, which leads to 



+ii 



V 



S_ ^5*4 



V 



< 1. 



Let us state what we have just proved into the following proposition. 



Proposition 10. The boundary condition (|2..Slj) is maximal non-negative if and only if 
there exists al- x Ij^ matrix M{x, t) such that 



E{x,t) 



( .T 



M{x,t) 



(2.46) 
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with 



+ll 



I0+I-+1 



M 



< 1. 



\ 



The following corollary is an immediate consequence of Proposition IIUI 

Corollary 11. A sufficient condition for having a maximal non-negative boundary condition 
of the form (|2.H1|) . with E{x,t) given by (j2.4(i|) . is 



2.4.2 Boundary Conditions for Constrained Systems 



In this subsection we establish some connections between the problems (|2.2()|) - (|2.22j) and 



(|2.23|) - (|2.24j) restricted to a bounded smooth domain 17. More precisely, we analyze the links 
between the boundary subspaces and boundary matrices associated to the two problems. 



First of all, it is easy to see that the boundary matrix associated to (|2.2,Sjl - (|2.24j> is given 
by 

' An{x) Bl{x) 

Bn{x) 

where, as in Subsection 12.4.11 An{x) = — X^j^i^j^"'; ^"^^ Bn{x) = —"^jLinjE^, with 
n(x) = (ni, . . . , riAr) the outer normal to 17 at x G dil. 



An{x) 



Lemma 12. 



for all x £ d^l. 



Anix){Ker: Bn{x)) C KerB„(x) 



Proof. Let v G Kei Bn{x) and u £ Ker B such that v = u{n). Applying the Fourier trans- 
form to (|2.9|) . it follows that Bn{x)An{x)u{n) = 0, or Bn{x)An{x)v = 0. So, ^„(a;)(Ker Bn{x)) C 
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KerB„(x), for all X G □ 
Lemma 13. 



Ker Anix) = (Ker A„(x) n Ker Bnix)) x Ker B'^{x), (2.47) 



for all X G dQ. 
Proof. Clearly, 



(Ker A„(x) n Ker Bnix)) x KerB^(x) C Ker A„(x), (2.48) 



for all X G dQ. 



Let I ^ I G Ker A„(x). Then 



Anu + z = 0, (2.49) 
BnU = 0. (2.50) 

From ()2.5U|) . u G Keri?„, and so, AnU G Keri?„, from the previous lemma. Applying B^ 
to (|2.49jl . it follows that B^B'^z = 0, which implies = 0, by multiplying it with . 

Thus z G KerB^. Returning to 1)2. 49(1 . observe that u G Kerj4„. Putting together all this 
information, we have that 



u 



T 



G (Ker An n Ker S„) x Ker 



Hence, 

Ker]4„(x) C (Ker A„(x) n Ker B„(x)) x Ker5^(x), (2.51) 

for all X G dO.. 

Prom (ITIHl) and (tTTni . we get (tTTTD . □ 
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Lemma 14. If u is a non-negative vector for An and z-LBnU, then 



u 



z 



) 



is non-negative 



for An. 



Proof. Under the given hypotheses, the conclusion is a simple consequence of the following 



An immediate corollary of this lemma is 

Corollary 15. The subspace N = N x [Bn{N)]-^ is non-negative for An if and only if N 
is non-negative for An. 

Inhomogeneous Boundary Conditions 

In this part, we consider the problem of finding u{x, t) G M"*, x G $7 C M^, t > 0, for l|2.2Uj) . 
subject to initial condition (|2.22|) . constraints (|2.2H) . and linear inhomogeneous boundary 
conditions 



Assume 5 is a vector function defined on 517 for all time t > such that there exists g 
satisfying the constraint equation (|2.21j) and g = Eg on the boundary for all time t > 0. 
Then, by substituting u = u — g, we arrive to the constrained initial homogeneous boundary 
value problem 



equalities 




□ 



E{x,t)u{x,t) = g{x,t) on dil x (0,r). 



(2.52) 



u = Au + A~g-~g + f, 



(2.53) 



n(x,0) = uq{x) -g, 



(2.54) 



Bu = 0, 



(2.55) 



Eu = 0. 



(2.56) 
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It is easy to see that the compatibihty conditions for this last problem are satisfied and, in 
fact, ()2.53|) - (|2.56|) is equivalent to the original constrained initial inhomogeneous boundary 
value problem (|2.20jl ~ (|2.22|) and (|2.52)) . Therefore, in this way and for a restricted set of 
boundary inhomogeneities, the treatment of inhomogeneous boundary conditions reduces 
to the treatment of the corresponding homogeneous ones. 
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Chapter 3 

Hyperbolic Formulations of 
Einstein Equations 

3.1 Introduction 

Einstein's equations can be viewed as equations for geometries, that is, its solutions are 
equivalent classes under spacetime diffeomorphisms of metric tensors. To break this diffeo- 
morphism invariance, Einstein's equations must be first transformed into a system having 
a well-posed Cauchy problem. The initial method to solve this problem has been by "fixing 
the gauge" |2S]) or in other words, by imposing some conditions on the metric com- 
ponents which would select only one representative from each equivalent class of Einstein's 
solutions. By using an ingenious choice of gauge fixing, the so called "harmonic gauge," 
Einstein's equations can be converted to a set of coupled wave equations, one for each 
metric component. Thus, by prescribing at an initial hypersurface values for the variables 
and their normal derivatives, one gets unique solutions to this system of wave equations. 
In Appendix B, we explain this approach, restricted to the linearized case for the sake of 
simplicity. For a more detailed discussion on this topic, see jlj, JH]) and references therein. 

Another way to deal with the gauge freedom of Einstein's equations is by prescribing the 
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time foliation along evolution, that is, by prescribing a lapse-shift pair along evolution 0. 
Einstein's equations are then decomposed into evolution equations and constraint equations 
on the foliation hypersurfaces. An analogous decomposition occurs for Maxwell's equations, 
which are canonically split into constraint (divergence) equations and evolution (curl) equa- 
tions. Both the constraints and the evolution equations are not uniquely determined by this 
procedure: by taking combinations of the constraints or/and by adding any combination of 
constraints to any of the evolution equations one gets an equivalent decomposition. 

There have been numerical schemes to solve Einstein's equations based on the ADM de- 
composition ini, but they have had only a very limited success. The instabilities observed 
in the ADM based schemes might be at least in part caused by the weakly hyperbolicity 
of the first order differential reduction of the ADM evolution equations, as argued in f36j. 
In fact, it is known that some of the stability problems of numerical schemes are due to 
properties of the equations themselves. By rewriting the equations in a different form, one 
can obtain better stability of computations for the very same numerical methods. 

There have been a large number of first order hyperbolic formulations derived from the ADM 
decomposition in recent years [^j, 0, [SSI, 0' 0) Ml; |HEj; IHji among others. To give a 
survey of all of them appears to be a very difficult and extensive task, which is out of the 
scope of this dissertation (see |28l for a more comprehensive review). All such formulations 
must be equivalent since they describe the same physical phenomenon. However, they can 
admit different kinds of unphysical solutions which can grow rapidly in time and overwhelm 
the physical solution in numerical computations. This is one reason (signaled in [HS]) why 
some formulations of Einstein's equations behave numerically better than others. From this 
point of view, first order symmetric hyperbolic formulations of Einstein's equations present 
a special attraction for a couple of reasons. First of all, there is a large body of experience 
and numerical codes that are stable for numerical simulations of FOSH systems derived 
from various applications (transport equations, wave equations, electromagnetism, etc.). 
Also, the symmetric hyperbolicity of the system ensures well-posedness and gives bounds 
on the solution growth. In fact, hyperbolicity refers to algebraic conditions on the principal 
part of the equations which imply well-posedness, that is, if appropriate initial data is given 
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on an appropriate hypersurface, then a unique solution can be found in a neighborhood of 
that hypersurface, and that solution depends continuously, with respect to an appropriate 
norm, on the initial data. 

In this chapter, we discuss some FOSH formulations of Einstein's equations introduced in 
recent years. Moreover, a new formulation due to Arnold 2, is presented and analyzed. 

3.2 Einstein Equations 

In general relativity, spacetime is a 4-dimensional manifold of events endowed with a pseudo- 
Riemannian metric 

ds^ = Qapdx^dx^ . (3.1) 

This metric determines curvature on the manifold, and Einstein's equations relate the cur- 
vature at a point of spacetime to the mass-energy there: 

Gap = SnTafs, (3.2) 

where Ga/3 is the Einstein tensor and T^/s is the energy-momentum tensor. 

The Einstein tensor G is a second order tensor built from the given metric as follows: 

Ga/B = Ral3 — -^RgafB- 

Here Ra/s is the Ricci tensor: 

Raf3 = RaSjB^ 

where 

is the Riemann curvature tensor. 
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The Christoffel symbols are defined as 
By R we denote the scalar curvature 

R = 

The energy-momentum tensor T^p can be better understood by looking at two of the sim- 
plest energy-momentum tensors in general relativity, namely, the energy-momentum tensors 
for incoherent matter or dust and for a perfect fluid. 

a) Incoherent matter (non-interacting incoherent matter or dust) 

Such a field may be characterized by two quantities, the 4-velocity vector field of flow 

dx"- 

where r is the proper time along the world-line of a dust particle and a scalar field 

Po = Po{x) 

describing the proper density of the flow, that is, the density which would be measured by 
an observer moving with the flow (a co- moving observer). 

The simplest second-rank tensor we can construct from these two quantities is 

and this turn out to be the energy-momentum tensor for the matter field. 

Now, let us investigate this tensor in special relativity in Minkowski coordinates. In this 
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case 



= — = 7(1, Ux,Uy, Uz), 

where 7 = (1 — n^)~^/^ and r is the proper time defined by 



dT^ = ds^ = dt^ - dx^ - dy^ - dz^ = dt\l - u^) = -f'^dt^ 



Then, the T°° component of T"" is 



ab 



-t — Po 



dx° dx^ dt^ 
Po- 



dr dr 



'dT^ 



This quantity has a simple physical interpretation. First of all, in special relativity, the mass 
of a body in motion is greater than its rest mass by a factor 7 (m = 77720). In addition, if we 
consider a moving three-dimensional volume element, then its volume decreases by a factor 
7 through the Lorentz contraction. Thus, from the point of view of a fixed as opposed to 
a comoving observer, the density increases by a factor 7^. Hence, if a field of material of 
proper density po flows past a fixed observer with velocity u, then the observer will measure 
a density p = 7^po- 

The component T^^ may therefore be interpreted as the relativistic energy density of the 
matter field since the only contribution to the energy of the field is from the motion of the 
matter. 

The components of T"** are 



/ 1 



U-r 



Uz 



2 

'^y '^x'^y '^y '^y'^z 



uz 



Next, we will show that the equations governing the force- free motion of a matter field of 
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dust can be written in the following very succinct way 



dbT"' = 0. 

When a = 0, this equation becomes exactly the classical equation of continuity 

dp 

— + dtv [pu) = 0, 

which expresses the conservation of matter with density p moving with velocity u. Since 
matter is the same as energy, it follows that the conservation of energy equation for dust is 

dbT'^'' = 0. 



Writing the equations corresponding to a G {1, 2, 3}, we get 



d , , d . , d , , d , 



Combining this with the equation of continuity, we obtain 

which is the Euler equation of motion for a perfect fluid in classical fluid dynamics (in the 
absence of pressure and external forces). 

We have seen that the requirement that the energy-momentum tensor has zero divergence 
in special relativity is equivalent to demanding conservation of energy and conservation of 
momentum in the matter field (hence the name energy-momentum tensor). 

If we use a non-flat metric, then the conservation law 

dhT'^'' = 
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is replaced by its covariant counterpart 

VbT"^ = 0. 

b) Perfect Fluid 

A perfect fluid is characterized by three quantities 

1. a 4--velocity u"' = dx°'/dT, 

2. a proper density po = po{x), 

3. a scalar pressure p = p{x). 

Observe that, if p vanishes, a perfect fluid reduces to incoherent matter. This suggests that 
we take the energy-momentum tensor for a perfect fluid to be of the form 

T"^ = pou^'u^ + pS"'^ 

for some symmetric tensor S""^. Since this tensor depends on the velocity and the metric, 
the simplest assumption we can make is 

where A and /x are constants. 
Considering the conservation law 

SbT'^* = 

in special relativity in Minkowski coordinates and demanding that it reduces in an appro- 
priate limit to the continuity equation and the Euler equation in the absence of body forces, 
we obtain that A = 1 and jj, = —1. Therefore, the energy-momentum tensor of a perfect 
fluid is 

T'^'' = {po+pKv!'-pg''\ 
34 



In the full theory, we again take the covariant form VbT = for the conservation law. In 
addition, p and p are related by an equation of state which, in general, is an equation of the 
form p = p{p, T), where T is the absolute temperature. Usually, T is constant and so, that 
equation of state reduces to p = p{p)- 

The Einstein equations can be viewed in three different ways: 

1. The field equations are differential equations for determining the metric tensor gab from 

a given energy-momentum tensor Tab- An important case of the equations is when Tab = 0, 
in which case we are concerned with finding vacuum solutions. 

2. The field equations are equations from which the energy-momentum tensor can be read 
off corresponding to a given metric tensor gab- In fact, this rarely turns out to be very 
useful in practice because the resulting Tab are usually physically unrealistic. In particular, 
it frequently turns out that the energy density goes negative in some region, which is rejected 
as unphysical. 

3. The field equations consist of ten equations connecting twenty quantities (the ten com- 
ponents of gab and the ten components of Tab)- I^i this way, the field equations are viewed 
as constraints on the simultaneous choice of gab and Tab- This point of view is useful 
when one can partly specify the geometry and the energy-momentum tensor from physical 
considerations and then the equations are used to determine both quantities completely. 

Remark 4. 1. Each of the 10 equations Gap = SttTq^ is a 2nd order PDE in 4 inde- 
pendent variables and 10 unknowns. These 10 equations involve lots of terms. 

2. The equations are not independent, due to the Bianchi identities VaG"^ = 0. The 
energy-momentum tensor also must satisfy these identities. 

3. Solutions are not unique, because of gauge freedom ( any diffeomorphism of the mani- 
fold gives a reparametrization, and hence another solution). 

A subtle consequence of Einstein's equations is that relatively accelerating bodies emit 
gravitational waves. These gravitational waves are very slight variations in the spacetime 
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Figure 3.1: A spacetime diagram for the ADM decomposition illustrating the definition of 
the lapse function, N, and shift vector 

metric tensor, which propagate at the speed of light. It is generally believed that extremely 
violent movements of huge masses, such as collisions of black holes, should generate de- 
tectable gravitational waves. However, because of their tiny amplitude, gravitational waves 
have eluded detection until now. When gravitational waves will be detected, we will must 
determine the cosmological event that could have caused them. This is in fact an inverse 
problem and, as usual, we need the solution of the direct problem, which in this case is the 
numerical solution of Einstein's equations. 



3.3 ADM (3+1) Decomposition 



In numerical relativity, the Einstein equations are usually solved as an initial-boundary 
value problem. In other words, the spacetime is foliated and each slice is characterized 
by its intrinsic geometry 'jij and extrinsic curvature Kij. Subsequent slices are connected 
via the lapse function and shift vector /3* (see Figure IXTjl . The ADM decomposition [H] 
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(also (53j) of the line element 

ds'^ = -N^dt^ + jij{dx' + (3'dt){dx^ + p^dt), i, j = 1, 2, 3, 

allows one to express six of the ten components of Einstein's equations in vacuum as a 
system of evolution equations for the metric jij and the extrinsic curvature Kij: 

% = -2NK,,+2V(^,(3j), (3.3) 

kij = N[R,, + {K\)Kij - 2KiiK^j] + (3^ViKij + KuVj(3^ + Ki,Vi(3^ - V,V,-7V, (3.4) 

m + {Klf - K,,K^^ = 0, (3.5) 

V^Kij - ViK] = 0, (3.6) 

where we use a dot to denote time differentiation. The spatial Ricci tensor R has components 
given by second order spatial partial differential operators applied to the spatial metric 
components ^ij. Indices are raised and traces taken with respect to the spatial metric, and 
parenthesized indices are used to denote the symmetric part of a tensor. 

The system of equations for and Kij is first order in time and second order in space. 
It is not hyperbolic in any usual sense, and direct numerical approaches have been unsuc- 
cessful. Therefore, many authors have considered reformulations of (|3.3|) . (|3.4|) into more 
standard first order hyperbolic systems. Typically, these approaches involve introducing 
other variables, like the first spatial derivatives of the spatial metric components ^ij (or 
quantities closely related to them) . In the rest of this chapter we present some of the most 
important first order hyperbolic formulations derived from the ADM decomposition, as well 
as their linearizations around the Minkowski spacetime. In the last section of this chapter 
we introduce a new first order symmetric hyperbolic formulation of the linearized ADM 
equations around Minkowski's spacetime that has surprising resemblances with Maxwell's 
equations. This motivates the introduction of the linearized ADM equations around the 
flat spacetime in the following subsection. 
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3.3.1 Linearized ADM Decomposition 



In this subsection, we derive the hnearized ADM decomposition around the Minkowski 
spacetime. 

A trivial solution to the ADM system (|3.,'-{j) " (|3.Hj) is the Minkowski spacetime in Cartesian 
coordinates, given by ^ij = 6ij, Kij = 0, /?* = 0, = 1. To derive the linearization of (|3.8j) - 
H3.6|) about this solution, we write jij = Sij + Qij, Kij = Kij, j3''' = 13^ , N = 1 + N , where the 
bars indicate perturbations, assumed to be small. If we substitute these expressions into 
H3.3|l - H3.6|l and ignore terms which are at least quadratic in the perturbations and their 
derivatives, then we obtain a linear system for the perturbations. Dropping the bars, the 
system is 



The usual approach to solving the system p.Tjl - p.lUf) is to begin with initial data gij{0) 
and Kij{0) defined on and satisfying the constraint equations ()3.9() . 1)3. 1U() . and to define 
Qij and Kij for t > via the Cauchy problem for the evolution equations (|3.7|) . (|3.8|) . It 
can be easily shown that the constraints are then satisfied for all times. Indeed, if we apply 
the Hamiltonian constraint operator defined in 1)3. 9() to the evolution equation ()3.7|) and 
apply the momentum constraint operator defined in (|3.1Ujl to the evolution equation 1)3. 8j) . 
we obtain the first order symmetric hyperbolic system 



gij = -2Kij + 2(9(i/?j), 
Kij = d^d(^jgi)i - ^d^digij - ]^didjg\ - didjN, 
C:=d^id'gij-d,gl) = 0, 



(3.7) 



(3.8) 



(3.9) 



Cj := &Kij - djKl = 0. 



(3.10) 




Thus if C and Cj vanish at t = 0, they vanish for all time. 
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3.4 Hyperbolic Formulations 



In this section, we present a number of popular hyperbolic formulations of Einstein's equa- 
tions derived from the ADM decomposition. In the last section, we present and analyze a 
new first order symmetric hyperbolic formulation of Einstein's equations in the linearized 
case due to Arnold |2j. 

3.4.1 Kidder-Scheel-Teukolsky (KST) Family 

In this subsection, we present a many-parameter family of hyperbolic representations of 
Einstein's equations intoduced by Kidder, Scheel, and Teukolsky [HS] . 

In order to write the evolution equations and ()3.4() in first-order form, we have to 

eliminate the second order derivatives of the spatial metric. For this purpose, we introduce 
new variables 

dkij = dkjij- (3.11) 

Since 

^ij = ^7"^(^(j9a7bj) + c?a5(j7j> - dadbjij - d(^idj)-fab) + lower order, 
the evolution system 1)3. 3|1 . 1)3. 4|) . together with ()3.11|) give 

^ij = -2NKij + ViPj + VjPi, (3.12) 
^ij = ^^7''*(%4fei) + dad^ij)b - dadbij - d(idj)ab) - didjN + (3.13) 
dk^J = -2NdkK,, - 2K,,dkN + dkiViPj + VjPi). (3.14) 

Since we have introduced a new variable that we will evolve independently of the metric, 
we get additional constraints, 

Ckij ■= dkij - dklij = 0, (3.15) 
Ckiij '■= d[kdi]ij = 0. (3.16) 
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The system ()3.12() - ()3.14() has been proven to be only weakly hyperbolic. Its characteristic 
matrix has eigenvalues {0, ±1}, but does not have a complete set of eigenvectors. Fortu- 
nately, the hyperbolicity of the system can be changed by densitizing the lapse and adding 
constraints to the evolution equations. 

We densitize the lapse by writing: N = j'^e^, where 7 = det {jij), a is the densitization 
parameter, which is an arbitrary constant, and Q is the lapse density, which will be chosen 
independent of the dynamical fields. 

By adding terms proportional to the constraints, we can modify the evolution system 1)3.12(1 - 
(|3.14|) without affecting the physical solution: 

k.,j = (...) + B^iijC + CA^7"'C'„(,,)fe, (3.17) 
dkij = (•••) + "n^lHiC^) + x^laCk, (3.18) 

where (. . .) represents the same thing that was before, C is the Hamiltonian constraint, Ci 
are the momentum constraints, and {0, 77, x} ^-re arbitrary parameters. 

By carrying out the computations, the evolution system, up to the principal part, is now 
given by (System 1 in |36j ) : 

iij = o> 

Kij = lNr'[dadki, - (1 + 0^ad(^^J)b " (1 " 0^(^dabj) + 

(1 + 2a)a(i(i,),, - 0^^J7'''^ad,db + Oj^jJ^'^dadbcd], (^'l^) 

dkij = -2NdkKij + 7V7"^(r/7fc(j9aKbj) + xHjdaKbk - V7k{idj)Kab- 

-XlijdkKab)- 

The eigenvalues of the characteristic matrix of this system are {0, ±1, =bci, =bc2, ±£3}, where 
ci = \/2^, C2 = 2-3/Vr/ - 4r/cj - 2x - 12o-x - 3r/C, and 
C3 = 2-1/V2 + A0 -r]-2er] + 2x + ^Bx - r]Q. 

If all Ci are real, then it can be proven that the system is strongly hyperbolic unless one of 
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the following conditions occurs: 



Ci = 0, (3.20) 
ci=C3^ 1, (3.21) 
ci = C3 = 1 7^ C2. (3.22) 

By differentiating the constraints in time, we obtain the following system for the evolution 
of the constraints (up to the principal part) 

C- = -^(2 - ry + 2x)NY'dpC,, (3.23) 



Ci = -(1 + 2e)NdiC+^NjP'iY'[i'i- - OdgCprsi + (1 + OdpCsigr - (1 + 2a)dpCgirs], (3.24) 



Ckij = 0, (3.25) 

Ckiij = ^v]^{ljiidk]Ci + Jiiidk]Cj) + xNjijd[kCiy (3.26) 

The eigenvalues for the characteristic matrix of this last system is a subset of the eigenval- 
ues of the evolution equations {0, ±02, ±03}. Moreover, the constraint evolution system is 
strongly hyperbolic whenever the regular evolution system is strongly hyperbolic. 

We define two new variables: the generalized extrinsic curvature and the generalized 
derivative of the metric M^ij 



Pij = Kij + zjijK, so Kij = Pij - Y:p^7ij-P, (3.27) 

^1 



^kij = hkdkij + ed^ij)k + lij{adk + hhk) + 7k(i{cdj) + dbj))], (3.28) 



where we introduce seven additional parameters {a, b, c, d, e, k, z}, and dj = -y djab, bj 

r'dabj. 
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The System 1, (|3.19() . rewritten for these new variables, up to the principal parts, becomes 
(System 2 in 

iij = 0, 

Pij = -Nj^'^ifiidaMbij + fi2daM^ij)i, + fi3d(^iM^f,j) + iJ-4d^iMj)^k+ 

+l^5lij7"'daM,db + f^67ij7"'daMh,d), (3-29) 
+^^l''^lk{idj)Pab + y^mjl^^dhPab), 

where jn = iii{a, 9, (, rj, x, a, b, c, d, e, k, z) and i^i = Vi{a, 9, (, rj, x, a, b, c, d, e, k, z). 

The System 2, (|3.29|) . is also strongly hyperbolic. Moreover, it has the same eigenvalues as 
System 1, (|3.19p . but the eigenvectors are different. 

Conclusion: By densitizing the lapse, adding constraints to the evolution equations, and 
changing variables, Kidder, Scheel, and Teukolsky [HH] got a twelve-parameter family of 
strongly hyperbolic formulations of Einstein's equations. They have observed that the choice 
of parameters can have a huge impact on the amount of time and accuracy of numerical 
simulations. Unfortunately, nobody knows why one particular parameter choice behaves 
much better than others. 

Some well-known formulations can be recovered by making appropriate choices for parame- 
ters. For example, we can recover the Fritteli-Reula (FR) system |25j if the following choice 
of parameters is taken in (|3.29p : 

{a, 9, C, rj, X, z, k, a, b, c, d, e} = {1/2, -1, -1, 4, -2, -1, 1, -1, 0, 0, 0, 0}. 

In the next subsection, we present another popular hyperbolic formulation which can be 
recovered from (|3.29|) . 
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3.4.2 Einstein-Christoffel (EC) formulation 

The EC formulation was originally derived directly from the ADM system (|c{.3|) - (|c{.6p by 
Anderson and York [7] in 1999. It can be also recovered if we make the following choice of 
parameters in (|3.29|) : 

{a, 9, C, ri, X, z, k, a, b, c, d, e} = {1/2, 0, -1, 4, 0, 0, 1, 0, 0, 2, -2, 0}. 

In this case, the system is (up to the principal parts): 



Essentially the coupled part of this system, i.e. the last two equations, is a set of six coupled 
quasilinear scalar wave equations with nonlinear source terms. 

Linearized Einstein ChristofTel Formulation 

First, we replace the lapse in (|3.3|) - (|3.6|) with where a denotes the lapse density. 
A trivial solution to this system is Minkowski spacetime in Cartesian coordinates, given by 
7jj = 5ij, Kij = 0, = 0, a = 1. To derive the linearization, we write 'yij = 6ij + gij, 
Kij = Kij, = a = 1 + a, where the bars indicate perturbations, assumed to be small. 
If we substitute these expressions into ((3ISl)-(|Sini) (with N = a^), and ignore terms which 
are at least quadratic in the perturbations and their derivatives, then we obtain a linear 



dtjij = 0, 



(3.30) 



dtMkij = - 



NdkPij. 
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system for the perturbations. Dropping the bars, the system is 



Qij = -2Kij + 2d(^iPj), 
kij = d'-d(^jgi)i - ^d^digij - didjg\ - didja 



(3.31) 



(3.32) 



C:=d^{d'gij-d,gl) = 0, 
Cj := d'Kij - d,K\ = 0, 



(3.33) 



(3.34) 



where we use a dot to denote time differentiation. 

Remark 5. For the linear system the effect of densitizing the lapse is to change the coef- 



— 1/2 instead of —1, and the derivation of the linearized EC formulation below would not be 
possible. 

The usual approach to solving the system (|3.31|) - ()3.34|) is to begin with initial data gij{0) 
and Kij{0) defined on M'^ and satisfying the constraint equations 1)3. 33(1 . (|3.34|) . and to define 
gij and Kij for t > via the Cauchy problem for the evolution equations (|3.31|) . ()3.32|) . It 
can be easily shown that the constraints are then satisfied for all times. Indeed, if we apply 
the Hamiltonian constraint operator defined in 1)3. 33(1 to the evolution equation ()3.31|) and 
apply the momentum constraint operator defined in 1)3. 34() to the evolution equation 1)3. 32() . 
we obtain the first order symmetric hyperbolic system 



Thus if C and Cj vanish at t = 0, they vanish for all time. 

The linearized EC formulation provides an alternate approach to obtaining a solution of 
(|3.31jl - (|3.34jl with the given initial data, based on solving a system with better hyperbolicity 
properties. If gij, Kij solve 1)3. 31() - (|3. 34(1 . define 



ficient of the term didj 



jg\ in ()3.32|) . Had we not densitized, the coefficient would have been 



C 



2d^Ci 




fkij = ^[dkOij - {d^gii - dig\)5jk - {d^gij - djg\ 



(3.35) 
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Then —d^ fkij coincides with the first three terms of the right-hand side of H3.32() . so 

Kij = -d^fkij - didja. (3.36) 

Differentiating (|3.35|) in time, substituting (|3.3H) . and using the constraint equation (|3.34j) . 
we obtain 

fkij = —dkKij + Lkij, (3.37) 

where 

Lkij = dkd^iPj) - d^d[il3i]6jk - d^d[iPj]5ik. (3.38) 

The evolution equations 1)3. 36(1 and 1)3. 37(1 for Kij and fkij, together with the evolution 
equation (|3.31|) for gij, form the linearized EC system. As initial data for this system we 
use the given initial values of gij and Kij and derive the initial values for f^ij from those of 
gij based on 1)3. 35|) : 

fkijiO) = \{dkgij{0) - [d'guiO] - diglmSjk - [d'gij{0) - d^glmSik}. (3.39) 

A purpose of this dissertation is to study the preservation of constraints by the linearized 
EC system and the closely related question of the equivalence of that system and the 
linearized ADM system. Our interest is in the case when the spatial domain is bounded 
and appropriate boundary conditions are imposed, but first we consider the result for the 
pure Cauchy problem in the remainder of this subsection. 

Suppose that Kij and f^ij satisfy the evolution equations (|3.36jl and (|3.37|) (which decouple 
from (|3.31|) ). If K^j satisfies the momentum constraint ()3.34() for all time, then from 1)3. 36(1 
we obtain a constraint which must be satisfied by f^ij- 

d\d'fui,-djfJ) = Q. (3.40) 

The following theorem shows that the pair of constraints 1)3. 34() . 1)3. 4U() is preserved by the 
linearized EC evolution. 
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Theorem 16. Let initial data Kij(0), fkijiO) be given satisfying the constraints (|3.34() and 
(|3.4U() . Then the unique solution of the evolution equations ()3.36() . ()3.37|) satisfy (|3.34() and 
(|3.40|) for all time. 

Proof. It is immediate from the evokition equations that each component Kij satisfies the 
inhomogeneous wave equation 

Kij = d'^dkKij - d'^Lkij - didjd. 

Applying the momentum constraint operator defined in H3.34[) . we see that each component 
Cj satisfies the homogeneous wave equation 

Cj = d'^dkCj. (3.41) 

Now Cj = at the initial time by assumption, so if we can show that Cj = at the initial 
time, we can conclude that Cj vanishes for all time. But, from H3.36|l and the definition of 

Cj, 

C, = -dHd'fHj-d,f,/), (3.42) 

which vanishes at the initial time by assumption. Thus we have shown Cj vanishes for all 
time, i.e., holds. In view of ((2311), (IXIU)) holds as weh. □ 

In view of this theorem it is straightforward to establish the key result that for given initial 
data satisfying the constraints, the unique solution of the linearized EC evolution equations 
satisfies the linearized ADM system, and so the linearized ADM system and the linearized 
EC system are equivalent. 

Theorem 17. Suppose that initial data gij{0) and Kij{0) are given satisfying the Hamilto- 
nian constraint (|3.33() and momentum constraint (|3.34() . respectively, and that initial data 
fkiji^) is defined by ()3.39|) . Then the unique solution of the linearized EC evolution equa- 
tions (|m?T|) . (pUr7|) satisfies the linearized ADM system (pOT|) . 
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Proof. First we show that the initial data fkij{0) defined in ()3.35|) satisfies the constraint 
H3.4U() . Applying the constraint operator in (|3.4U|) to H3.35() we find 

which vanishes at time by H3.33() . From Theorem ll6l we conclude that Cj = for all time, 
i.e., (TOIl) holds. Then from (|OT|) and ^^HM we see that C = -2d^Cj = 0, and, since C 
vanishes at initial time by assumption, C vanishes for all time, i.e., (|3.33jl holds as well. 

It remains to verify ()3.32|) . From 1)3. 37(1 and ()3.31|) we have 

fkij = ijdkgij - d^d[iPii6jk - d^d]^il3j^6ik. 
Applying the momentum constraint operator to (|3.31|) and using 1)3.34(1 . it follows that 

so fkij - [dkgij - {d'-gii - dig\)5kj - {d'-gij - djg\)5ki]/2 does not depend on time. In view of 
(POU) . we have (lO^ . 

Substituting (|3.35|) in (|3.36|) gives 1(3. 32(1 . as desired. □ 
3.4.3 Alekseenko— Arnold (AA) Formulation 

In this subsection we present a first order symmetric hyperbolic formulation of the full 
nonlinear ADM system (|3.3() - (|3.6() which involves fewer unknowns than other hyperbolic 
formulations and does not require any arbitrary parameters. The hyperbolic systems in- 
volves 14 unknowns, namely the components of the extrinsic curvature and eight particular 
combinations of the first derivatives of the spatial metric. In order to derive the AA formu- 
lation, we introduce the notations S for the 6-dimensional space of symmetric matrices and 
T for the 8-dimensional space of triply-indexed arrays (wijk) which are skew-symmetric in 
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the first two indices and satisfy the cyclic property 



Wijk + Wjki + Wkij = 0. (3.43) 



Define the operators M : C°°{R^,8) C°°(M^M3), L : C°°(R3,S) ^ C°°(R3,T), and 
L* : C°°(M3,T) ^ C°°(R3,S) by 

(Mn), = 2Y'id[j,u,]g, (3.44) 

{Lu)ijk = d[iUj]k, (3.45) 

iL*v)i, = -7qnrjdpvP^'^'\ (3.46) 

respectively. Observe that the operators L and L* are formal adjoints to each other with 
respect to the scalar products {u,w) := f UijW^^dx and {v,z) := f Vij^z^^^dx on the spaces 
C"^(M^S) and C°°(M^T) respectively. 



Next, we introduce new variables 

hjk = -^[{L^)ijk + (M7)[,7,]fc]. (3.47) 

Observe that {fijk) belongs to T and so, it has eight independent components. 

Now we are able to write down the new formulation derived in j^I for the ADM system 

^doKi, = -N{L*f)ij + ... = N^^anAf^""""^ + • • • > (3.48) 



V2 
1 
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dofijk = [LiNK)]ijk + ■■■ = dyi{NK)j^k + ■■■■ 



Here Oq := dt — P^di is the convective derivative and the omitted terms are algebraic ex- 
pressions involving 7jj, their spatial derivatives d^^ij, Kij, the lapse A^, and the shift /3. 



48 



Linearized AA Formulation 



The linearized AA formulation provides an alternate approach to obtaining a solution of 
p.Tjl - p.lUf) with the given initial data, based on solving a system with better hyperbolicity 
properties. If Qij := 7ij and Kij := Kij solve (|c{.7() - ()3.1U() . define 

Afey = --^[5[fc7i]i + {MJ)[k^^]J], (3.49) 

where {Mj)i = d^-fu - da\. 

Then, proceeding as in we obtain an evolution system for Kij and X^ij 

-^Kij = d^Xk(ij) - -^didja, (3.50) 

— ^Afcjj = d^i^Ki-^j + r]kij, (3.51) 

where 

Vkij = -^{djdikPi] + d'^d^^(3k]6ij - d'^d^^(3i^6kj). (3.52) 

The evolution equations (|3.5U|) and (|3.51|) for nij and Xkij form the linearized AA system. 
As initial data for this system we use the given initial values of and Kij and derive the 
initial values for Xkij from those of ^ij based on 1)3. 49(1 : 

Afc,,(0) = --L[5[fc7(0)i],- + (M7(0))[fc5,],]. (3.53) 

As shown in , if 7 and k satisfy the ADM system and A is defined by 1)3. 49() , then k and A 
satisfy the symmetric hyperbolic system ()3.5U|) . 1)3. 51(1 . Conversely, to recover the solution 
to the ADM system from (|3.50() . 1)3. 51() . the same k should be taken at time 0, and A should 
be given by (|3.53|) . Having k determined, the metric perturbation 7 is defined as follows 
from (|^ 

1^3 = 7ii(0) - 2 f {Kij - %/?,)) (s) ds. (3.54) 
JO 
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3.4.4 Arnold (A) Formulation 



In this subsection, we present and analyze the hnearized case of a new first order symmetric 
hyperbohc formulation of Einstein's equations due to Arnold |2j. 

Notations and Identities 



We begin by introducing a number of notations and listing a couple of useful identities. Let 
7 be a symmetric 3 x 3-matrix function and v a vector field (we use underbars for 3- vectors, 
double underbars for 3x3 matrices). 



skw7 = (7 — 7"'")/2 
/ 



-V3 V2 

V3 —vi 

\-V2 Vi J 



Skwv 



ciirlv= idu2/dz — dus/dy, du^/dx — dui/dz, dui/dy — du2/da 



curli- 7 is the matrix whose rows are the curls of the rows of 7 



curie 7 is the matrix whose columns are the curls of the columns of 7 



curls 7 = (curlr 7 + curie 7)/2 



\7 V is the matrix whose rows are the gradients of the entries of v 



div 7 is the vector whose components are the divergences of the rows of 7. 



ev= [V V + (y v)^]/2 



i?7 = ediv7 - 1/2 A7 - 1/2V Vtr7 



M7 = div 7 — V tr 7 



M*v = — €v + (div v)5 {S being the 3x3 unit matrix) 
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For two vector fields a and b, denote by a 6 the following symmetric matrix function 
aQb = {ab^ + ba^)/2. 

The proofs of the following identities imply only elementary computations and we leave 
them to the reader: 

curlr* = curlr (3.55) 



curie* = curie 



curls* = curls 



(3.56) 



(3.57) 



curie 7 = (curlj. 7)^ (therefore curlg 7 is symmetric) (3.58) 

trcurlr7 = (3.59) 

Mcurlr7 = (3.60) 

M curls 7 = 1/2 curl M7 = -^MSkw(M7) (3.61) 

curls 2 = curlr 7 - ^ Skw(M2) = curie 7 + ^ Skw(M7) (3.62) 

skw curlr 2 = ~ ^ Skw(M7) (3.63) 

^7 = ^ curie curlr 7 + ^ (div Mj)d (3.64) 

divMeu = (3.65) 
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M V V = for any function N (3.66) 

MM*v = -divev-Vdivv (3.67) 
A New FOSH Formulation of Linearized ADM 

If we choose lapse A^, shift and hnearize the ADM system about the Minkowski's flat 
space, the perturbations 7, k of the metric and the extrinsic curvature, respectively, are 
symmetric matrix fields satisfying 

2 = -2K + 2e^, (3.68) 

K = Rj-YY^^ (3.69) 
divM7 = 0, (3.70) 
Mk = 0. (3.71) 

Of course, (|3.68|) - (|3.71jl is exactly the linearized ADM system ()3.7|) - ()3.1U|) written in vector- 
matrix notation. 

This system should be supplemented with initial conditions 

2(0) = 70, k(0) = Ko, (3.72) 

(that satisfy the constraints (|3.7fl|) and (|3.71|) ) and boundary conditions if the domain has 
frontier. 



Lemma 18. 



= V(divM7). (3.73) 



Proof. It follows from (|3.69)) . (|3.7fl)) . and the identity: 

div curie 7 = curl div 7. 
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□ 

Theorem 19. If the initial conditions 70 and kq satisfy the Hamiltonian constraint 1)3 .701) 
and the momentum constraint (|3.71|) respectively, then the constraints are automatically 
satisfied for all time by any solution of the pure Cauchy problem (|3.(i8j) . (|3.69|1 . and (|3.72j) . 

Proof. From (TTO . (TUiHl) . and (TTtHI) . observe that 

^^^(divM^) = -2divMK = A(divM2). 

Denote by (/> = divMQ;. Then = A(?i, 0(0) = 0, and 0(0) = 0. Therefore, 

= div M7 = 0. 

Moreover, from Lemma [THl 

Mk = V(divM7) = 0, 
and from Mkq = 0, we get Mk = for ah time. □ 

From (|I;BH1) and (|T7T|) . it fohows that 

d'^K = -2Rk - V V iV, (3.74) 
Mk = 0. (3.75) 

Taking into account ()3.62|) . ()3.64|) . and 1)3. 70(1 . the equation ()3.74|) transforms into 

dfK = - curls curls k - V V iV. (3.76) 

Introduce = and = curls k. Then the equations (|3.76|) . 1)3. 75p induce the following 
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first order symmetric hyperbolic system with constraints 



£ = -curls^- VViV, (3.77) 

/i = curls v, (3.78) 

= 0, (3.79) 

= 0. (3.80) 

Here, the initial data is 

1/(0) = ^7(0) - V V N{0), /x(0) = curls k(0). (3.81) 



If 7o and kq satisfy the Hamiltonian constraint (|3.7U() and the momentum constraint 1)3. 71(1 
respectively, then it is not hard to see that the compatibility conditions Mi/(0) = and 
Mfj,{0) = are satisfied. 

Proposition 20. Suppose that 70 and kq satisfy the Hamiltonian constraint (|3.7U|1 and the 
momentum constraint (|3.71() respectively. Then, the problems ()3.68|) - ()3.72|1 and (|3.77|) - 
1)3. 81() are equivalent. 

Proof. From construction, it is obvious that once we have a solution for (|3.68|) - (|3.72|) . then 
(TrrTll - iprKT]) has solution. 

Let us prove that the converse is also valid. Assume we can solve (|3.77|) - H3.81|) and define 

K = Ko + / u{s)ds, (3.82) 
~ ~ Jo ~ 

7 = 70 + / [-2k{s) + 2e(3{s)] ds. (3.83) 
~ ~ Jo 

Obviously, 7(0) = 70 and k(0) = kq and so, (|3.72|) is verified. Now, we prove that 1)3. 7U() is 
satisfied. From l|3.61|) it follows that 



dt{Mu) = --curlM^, 
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(3.84) 



and from here 

dtidivM^) = 0. (3.85) 
Since divMi/g = 0, we get divMu = for all time. 

Next, we show that k satisfies dSZZU). From ((SZZHl), (jMU), and (|TH^ it follows that 

/i = curls K (3.86) 

for all time. By using H3.77() - (|3.78() . (|3.82|) . and (|3.86l) we can see that Mk satisfies the 
initial value problem 

df{MK) = ^AMk, (3.87) 

Mk(0) = 0, dt{M^){0) = 0, (3.88) 

whose only solution is the trivial solution. So, Mk = for all time. 

Finally, let us prove that (|3.68|) and (|3.69|) are verified. Observe that (|3.68|) is trivially 
satisfied by taking the derivative of (|3.83|) with respect to time. The identity (|3.69)) is a 
little bit more dehcate, but it follows from (jSSll), (EZZZl), (EH, and (jiU^ 

= dtK = - curls curls k - VViV = -2Rk - VViV = dtiRj - VVA^), (3.89) 

and from (TUni . we get ^W^. □ 

Equivalent Unconstrained Initial Value Problem 

Observe that the constrained problem H3.77|) ~ (|3.81|) can be put into the form H2.10|) - (|2.11() 
for 

-curls \ / M \ / -VViV \ / uiQ) 

A=\ = M , i? = - , / = = - , and no = ' 

curls 0/ \0M/ \ 0/ \ /x(0) 
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Moreover, the condition H2.9|l is satisfied since 



Bx = 4^ Mfj, = Mu_ = 44> curls /i = curl^ ^ and curlg v_ = curlr v_ 4^ 
M curls ^ = M curls i/ = 44> BAx = 0. 



Therefore, we can transform our system into a 18 x 18 symmetric hyperbolic (unconstrained) 
system 
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The initial data for the above system reads 



P 

V 9 y 



(0) 



/ Rjo-VVNiO) 
curls kq 


\ 



(3.91) 
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Chapter 4 



Boundary Conditions for Einstein's 
Equations 

4.1 Introduction 

A very difficult situation is encountered if initial boundary value problems are considered 
for Einstein's equations. It has become clear in the numerical relativity community that in 
order for constraints to be preserved during evolution, the boundary conditions have to be 
chosen in an appropriate way. This problem is widely open in its full generality for Einstein's 
equations. The main difficulty comes from the necessity to deal with unphysical instabilities, 
which will always be triggered in numerical simulations by round-off errors and lead to 
constraint violating numerical solutions. In other words, the numerical solution of the initial 
boundary value problem with constrained initial data fails to solve the constraints soon after 
the initial time. This serious drawback motivates the quest for boundary conditions that 
preserve the constraints. Only in the last few years work along this line has been pursued; 
see , |14j . |15j . \2'A\ . j24j . [M\ . among others. In view of our work, of particular interest 
is the recent paper by Calabrese et al. J15' on the generalized [321 Einstein-Christoffel 
formulation [7j linearized around Minkowski spacetime, which employs techniques with 
some points in common with the ones used in our approach of the problem in this chapter. 
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It is here where the differential equations satisfied by the constraints are very important. In 
particular, if they form a hyperbolic system, then the study of its well-posedness is proven 
to be a very good starting point for what boundary conditions we must force upon the 
evolution system for the dynamical variables. By exploiting this idea, we have been able 
to develop a technique that provides maximal nonnegative constraint preserving boundary 
conditions for various hyperbolic formulations of Einstein's equations in the linearized case. 
This entire chapter represents our original contribution to the subject. 



In this section we consider a constrained first order symmetric hyperbolic system, which, as 
we shall prove, admits well-posed constraint preserving boundary conditions. The analysis of 
this model problem gives a good deal of insight for finding same type of boundary conditions 
for various hyperbolic formulations of Einstein's equations. Although the problem is simple, 
the techniques used here reveal our basic strategy to tackle the more complex case of 
Einstein's equations. 

For t > and x in M^, we are interested in finding a solution {wi^ Vi, Uij) for the following 
first order symmetric hyperbolic system 



4.2 Model Problem 



Wi = Vi, Vi = d^Uij + fi 




with the initial data 



(4.2) 



and subject to the constraint 



c ■= d'vi = 0. 



(4.3) 



Assume that the initial conditions (|4.2j) are compatible with the constraint (|4.3|) 



d'v^i = 0, = 0. 



(4.4) 
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Also, assume that the forcing terms satisfy the compatibihty condition 

d^fi = 0, for ah time t > 0. (4.5) 

Then, for the pure Cauchy problem (|4.1j) . (|4.2|) . it can be easily shown that the constraint 
H4.3|) is satisfied for all time. Indeed, if we differentiate twice in time the constraint ()4.3() 
and use the main system ()4.ip . it follows that 

c = d'v, = dt[d\d^uij + /,)] = d'd^iiij + dtid'fi) = d^djd'vi + dt{d'f,) = d^d,c + dt{d'u). 

(4.6) 

Since d^fi = for all time, it follows that C satisfies the wave equation on x M^. 

C = AC. (4.7) 
Moreover, from the compatibility conditions ()4.4|) 



So, 
Also, 



d'vi{0) = d'v^ = 0. 

C(0) = 0. 
d'vi{0) = d^d%{0)+d'M0). 



(4.8) 

(4.9) 
(4.10) 



The first term on the right side vanishes from the second compatibility condition in 1)4.4(1 . 
The second term vanishes because of (|4.5() . Thus, 

C{0) = 0. (4.11) 

From and (|irTT|) . we conclude that C = for all time. 
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Figure 4.1: A polyhedral domain. 



4.2.1 Constraint-Preserving Boundary Conditions for the Model Prob- 
lem 

In this subsection, we are interested in finding suitable boundary conditions for 1)4. 1() on 
the boundary dO, of a bounded spatial domain 17 such that the solution of the resulting 
initial- boundary value problem satisfies the constraint (|4.3jl for all time. 

On Polyhedral Domains 

First, we investigate the existence of constraint-preserving boundary conditions on polyhe- 
dral domains, with the result that there exists a set of such boundary conditions. Moreover, 
these boundary conditions are also maximal nonnegative, and so, the corresponding initial- 
boundary value problem is well-posed. 

Let be a polyhedral domain in M'^. Consider an arbitrary face of dQ and let n* denote its 
exterior unit normal. Denote by and /* two additional vectors which together n* form 
an orthonormal basis (see Figure ETT|) . The projection operator orthogonal to is then 
given by t- := niim^ + (and does not depend on the particular choice of these tangential 
vectors). Note that 

5l=n^+Tl TiT^ = Tt (4.12) 
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Consequently, 

ailj- = n^ajriiU + r^ajT-b^ for all vectors a;, (4.13) 
We shall prove that the following set of boundary conditions 

n'-n^Uij = 0, m^'Vi = 0, Vvi = 0, (4.14) 

is constraint-preserving and, together with (|4.1j) and (|4.2|) . leads to a well-posed initial- 
boundary value problem. Observe that these boundary conditions can be written as well: 

n'n^Uij = 0, T'jVi = 0. (4.15) 
Therefore, they do not depend on the particular choice of the tangential vectors m* and Z*. 

Theorem 21. Given initial conditions Wi{0), fj(0), ■Ujj(O) and forcing terms fi satisfying 
the compatibility conditions (|4.4() and ()4.5|) respectively, define Wi, Vi, and Uij for positive 
time by the evolution equations (|4.1() and the boundary conditions ()4.14|) . or 1)4. 15|) . Then, 
the constraint (|4.3j) is satisfied for all time. 

Proof. First we prove that C = on any boundary face of 0,. From the first identity of 
(|4.12|) . the constraint C can be decomposed as follows 

C = d'vi = n'n^djVi + T'^djVi. (4.16) 

By the third equation in 1)4. 1|) and the second identity of 1)4. 12() . 1)4. 16() reads 

C = n'n^Uij + T'^Ti^jV^. (4.17) 

From the boundary conditions 1)4. 15(1 . we know that n^n^Uij = for all time, and so the first 
term on the right-hand side vanishes. Similarly, we know that t^'^vi = on the boundary 
face, and so the second term vanishes as well (since the differential operator T^dj is purely 
tangential). We have established that C = holds on any boundary face of 

Finally, since C also satisfies the wave equation ()4.7)) . together with the trivial initial con- 
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ditions H4.9|l and H4.11() . we conclude that the constraint C vanishes for all time. □ 

Theorem 22. The differential system (|4.1|) . together with initial data (|4.2() and boundary 
conditions (|4.15() is well-posed. Furthermore, if the initial conditions Wi{<S), fj(0), Uij{Q) 
and the forcing terms fi are given satisfying the compatibility conditions (|4.4|) and (|4.5j) 
respectively, then the solution satisfies the constraint ()4.3() for all time. 

Proof. We begin by showing that the boundary conditions H4.15() are maximal nonnegative 
for the hyperbolic system (|4.H) . and so, according to the classical theory of 18 and [HZj 
(see also ^01), the initial-boundary value problem is well-posed. We recall the definition 
of maximal nonnegative boundary conditions. Let V := M'^ x x M^^^. Obviously, 
dimy = 15. The boundary operator An associated to the hyperbolic system (|4.1|) is the 
symmetric linear operator V — > V which assigns to any element {wi, Vi, Uij ) of V the 
element {wi, Vi, Uij) defined by 

Wi = 0, Vi = —n^Uij, Uij = —UjVi. (4.18) 

A subspace of y is called nonnegative for An if 

w'wi + v'vi + u'Hij > 0, (4.19) 

whenever (wj, Uij) € N and (t(7j, Wj, Uij) defined by (|4.18|) . We claim that the subspace 
N defined by the boundary conditions 1)4. 15(1 is maximal nonnegative. The dimension of 
N is clearly 12. S ince An has three positive, nine zero, and three negative eigenvalues, a 
nonnegative subspace is maximal nonnegative if and only if it has dimension 12, and so, the 
proof reduces to the verification of nonnegativity condition (|4.19|) . Since 

w^Wi + v^Vi + u^-'uij = —2v^n-'uij, (4.20) 

the verification of l|4.19jl reduces to showing that v^n^Uij < whenever (|4.15|) holds. In fact, 
we prove that v^n^Uij = if the boundary conditions 1)4. 15() are verified. For this, we use 
the orthogonal decomposition of = n^VkU^ + T^^Vk. Prom the second boundary condition 
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n' 




Figure 4.2: A 3D domain. 

in 1)4. 15() . the tangential part T^^Vk of vanishes. Thus, v^n^Uij = li^n^UijU^Vk- From the 
first boundary condition in (|4.15|) . the right-hand side vanishes, and so, v^n^Uij = 0. This 
ends the proof of the maximal nonnegativity of the boundary conditions (|4.15|) . 

The fact that the solution satisfies the constraint for all time follows from Theorem 12 II □ 
On Piecewise Regular Domains 

Here we investigate the existence of constraint-preserving boundary conditions on more gen- 
eral domains, namely on piecewise regular domains. We shall provide a set of such boundary 
conditions which generalize the boundary conditions (|4.15|) . Although these boundary con- 
ditions are not maximal nonnegative in general, and so, the classical theory of 18.: and 
j37j can not be employed, we have been able to prove an energy inequality which is a key 
ingredient in proving well-posedness. 

Let be a bounded domain in R^. Assume that the boundary of $7, denoted by dCl is 
an almost everywhere regular surface in M^, that is, for almost every point p £ dil, there 
exists an open set V in and a difeomorphism <p : U ^ V H 5^1, (x^, x^) — > (t){x^, x^), 
of an open set U C onto V H dQ C M^. The mapping (p is called a parameterization 
or a system of (local) coordinates in a neighborhood of p. From "HH^, Corollary 2, p. 183, 
we know that for all p £ there exists a parameterization in a neighborhood V p 
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such that the coordinate curves intersect orthogonahy for each q £ V (such a (p is called 
orthogonal parameterization). Since working in an orthogonal parameterization turns out 
to simplify our computations, we assume that (j) is an orthogonal parameterization. At each 
point p of V n do, define the tangent vectors ti = d(f)/dx'^, i = 1, 2, and the unit normal to 
the boundary (see Figure W?^ 

n= - -, . 4.21 

- IliAtal 

For the given parameterization (f>, the derivatives of n belong to the tangential plane and 
are given by 

dn 

(4.22) 

an 

-—=^ = 012*1 + «22i2- 

The trace of the matrix (ajj)i<j,j<2 does not depend on the choice of parameterization. 
Hence, the mean curvature H, which is one half of the trace of (ajj)i<j,j<2i is parameteri- 
zation independent. Since the principal curvatures ki and k2 are the roots of the quadratic 
equation 

an - k ai2 
det I I = 0, 

^21 022 — k 



we can write 



H = ^{au+a22) = ^iki+k2). (4.23) 



We shall prove that the following set of boundary conditions 

n'n^Uij + 2Hn'wi = 0, r'^Vi = 0, (4.24) 

is constraint-preserving and, together with (|4.1I) and (|4.2() . leads to well-posedness. 
First we prove the equivalent of Theorem 1211 under the new circumstances. 

Theorem 23. Given initial conditions Wi{0), fj(0), Uij{0) and forcing terms fi satisfying 
the compatibility conditions 1)4. 4|) and 1)4. 5|) respectively, define Wi, Vi, and Uij for posi- 
tive time by the evolution equations ()4.1)) and the boundary conditions ()4.24|) . Then, the 
constraint ()4.3)) is satisfied for all time. 
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Proof. Let E = {ti, ti), F = {ti, t_2) = 0, and G = {t_2, t_2) be the coefficients of the 
first fundamental form in the basis { ti, t2 }■ The derivatives of the vectors ti and t2 in the 
basis { ti, t2, n} are (see p. 232): 

^=T\,t_,+Ti,t_2 + L,n, 
ll = & = ryi+r?,t, + L,n, (4.25) 



dx"^ dx 

^22*1 + ^22*2 + L^n, 



dt_2 _ 

^-i22J:1^^22i 



where T^^ are the Christoffel symbols and Li = —auE, L2 = —ai2E = —a2iG, and 
L3 = —a22G. 

Since the constraint C satisfies the wave equation (j4.7|) with zero initial conditions (|4.9j) . 
(|4.11() . the proof of this theorem reduces to prove that C = a.e. on the boundary dfl for 
all time. From the first identity of (|4.12jl . observe that on the boundary the constraint C 
is d^Vi = d^^djVi = n^n^djVi + r^^djVi. From the main system (j4.1|) . the first term of the 
right-hand side is equal to n^n^iiij. Since the parameterization is orthogonal, the tangential 
projection of djVi is 

where t^, 1 < j < 3, represent the components of the tangent vector tj, z = 1,2 . By the 
Leibnitz rule of differentiation, the right-hand side of this identity can be written as follows 

r'^djv^ = ^t{dj{Avi) + ^id^^Vi) - ^Vit{djA - y4djti. (4.26) 

The first two terms of the right-hand side vanish due to the second boundary condition in 
(|4.24() combined with the tangentiality of the differential operators tjdj, i = 1,2. Further- 
more, by the chain rule, t-'^,djt\. = dt\/dx^ , k = 1,2. Thus, 



1 rlf^ 1 rlf^ 
r'^d,v. = -lv.^-U^,. (4.27) 
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By the first identity in ()4.25l) . the first term of the right-hand side of (|4.27|) is 

1 dt^ 1 1 

-E'^'dx^ = -E^iii^'^'^i) - E^uii'^^^) + (4-28) 

The first two terms of the right side of the last identity vanish due to the second boundary 
condition in (|4.24l) . and so 

1 dti 

= «ii^V. (4.29) 



Similarly, 



1 df 

-^Vi^ = a22n'vi. (4.30) 



Thus, 

r'^djVi = (an + a22)n'vi. (4.31) 

Therefore, d'^Vi = n^n^ iiij + {aii + a22)n'^Vi. By the first equation in (|4.1j) . Wi = Vi, combined 
with the definition of the mean curvature, H = (aii+a22)/2, it follows that d'^Vi = (n^n^Uij + 
2Hn^Wi)'. Finally, from the first boundary condition in 1)4.24(1 . C = d^Vi = a.e. on dQ for 
all time. This ends the proof of this theorem. □ 

Theorem 24. Let Q be a bounded domain in with (a.e.) regular boundary and the 
mean curvature H > almost everywhere. Given the initial conditions w^, vf , u^j and 
the forcing terms fi satisfying the compatibility conditions (|4.4() and ()4.5() respectively, the 
system of differential equations ()4.1|) , together with the boundary conditions (|4.24|) and the 
initial conditions ()4.2() . is well-posed. Moreover, its solution satisfies the constraint ()4.3() 
for all time. 

Proof. For any solution Wi, Vi, Uij of (|4.1I) satisfying the boundary conditions (|4.24() . define 
the energy 

E{t) = -[ I {wiw' + Viv' + Uiju'^) dx + 2 H{win^f ds\. 
2 Jn Jan 

Differentiating in time, we obtain 



E{t) = / {wiw'' + ViV^ + iiiju''^) dx + 2 I H{win'){win'') ds. 
Jn Jan 
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By the system (|4.1() . the derivatives in time can be substituted for spatial derivatives, and 
so 

E{t) = [ [viw' + (d^Uijy + Uv' + {djVi)u'^ dx + 2 [ H{v,n'){win') ds. 
Jn Jan 

Integrating by parts, it follows that 

E{t) = / {viw' + fy) dx + / v'n^Uij ds + 2 / H{vin''){win') ds. 
Jn Jan Jan 

By the orthogonal decomposition of = n'^n^vi + t^^vi, we obtain v'^n^Uij = n^n^UijU^vi + 
r'^^Vin^Uij. In view of the second boundary condition in (|4.24|1 . the second term of the 
right-hand side vanishes. Thus, v^n^Uij = n^n^UijuI'vi, and so 

E(t) = / {viw'' + fy) dx+ I {n'n^Uij + 2Hwin'){n'-vi) ds. 
Jn Jan 

Due to the second boundary condition in (|4.24jl . the boundary integral vanishes. Therefore, 

E{t) = [ {viw' + fy)dx. 
Jn 

Now, since H > (a.e.), it is easy to see that 

E{t) < E{t)+e{t), 

where 9{t) = {J^ /j/*dx)/2 for all positive time t > 0. By the Gronwall inequality, E{t) < 
e^[E{0) + Jq e~'^9{s) ds]. This energy estimate is a key ingredient in proving well-posedness 
by Galerkin approximations. 

The fact that the solution satisfies the constraint (|4.3|) for all time follows from Theorem 1231 

□ 
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4.3 Einstein-Christoffel Formulation 



4.3.1 Maximal Nonnegative Constraint-Preserving Boundary Conditions 

In this section, we provide maximal nonnegative boundary conditions for the hnearized EC 
system which are constraint-preserving in the sense that the analogue of Theorem^] is true 
for the initial-boundary value problem. This will then imply the analogue of Theorem 1171 
We assume that is a polyhedral domain (see Figure IXT|) and use the notations introduced 
in Subsection 14.2.11 

First we consider the following boundary conditions on the face: 

n'm^Kij = n'VKij = n^n'n^ hij = n^m'm^ hij = n'^l'l' f^j = n^'m'P hij = 0. (4.32) 

These can be written as well: 

nW^''K,j = 0, n^^n'ri^ h^J = 0, n W^/^i,- = 0, (4.33) 

and so do not depend on the choice of basis for the tangent space. We begin by showing 
that these boundary conditions are maximal nonnegative for the hyperbolic system (|3.36j) . 
(|3.37|) . and so, according to the classical theory of 18^ and |37j . the initial-boundary value 
problem is well-posed. 

Let V denote the vector space of pairs of constant tensors (Kij, fkij) both symmetric with 
respect to the indices i and j. Thus diml/ = 24. The boundary operator An associated to 
the evolution equations (|3.3(ij) . (|3.37|) is the symmetric linear operator V ^ V given by 

Kij = n^fkij, fkij = rikKij. (4.34) 

A subspace of y is nonnegative for An if 

KijK'' + fkijf^'' > (4.35) 
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whenever {Kij,fkij) G N and {Kij,fkij) is defined by (|4.34() . The subspace is maximal 
nonnegative if also no larger subspace has this property. The eingenvalues A of An, together 
with the corresponding eigenvectors {Kij, fkij), are the following: 

A = with multiplicity 12 and eingenvectors: (0, nif^ninj), (0, rrii^ni^imj'j), (0, m^n^j/j-)), 
(0, mkmiiTij), (0, nikm^^ilj-j), (0, nikklj), (0, IkniUj), (0, Ikn^nij-)), (0, lkn(^ilj)), (0, kminij), 
(0, km^ilj-j), (0, khlj), 

A = —1 with multiplicity six and eingenvectors: {niUj, —nkniTij), {rriimj, —nkmirrij), 
{klj, -rikhlj), {n^inij-j, -nkn(^mj)), {n(^lj), -nkU^ilj)), {l(imj), -nkmirrij)), 

A = 1 with multiplicity six and eingenvectors: {riinj, ni^niUj), {rriirrij, nkmirrij), {Ulj, nkhlj), 

Since An has six positive, 12 zero, and six negative eigenvalues, a nonnegative subspace is 
maximal nonnegative if and only if it has dimension 18. Our claim is that the subspace 
defined by H4.32() is maximal nonnegative. The dimension is clearly 18. In view of (|4.34() . 
the verification of (|4.35|) reduces to showing that n^fkijK^^ > whenever (|4.32() holds. In 
fact, n^fkijK^-^ = 0, that is, fkij and Kij are orthogonal (when 1)4.32(1 holds). To see this, 
we use orthogonal expansions of each based on the normal and tangential components: 

Kij = n^nin"'njKim + n^niT^Kim + T^n'^njKira + T^T^Kim, (4.36) 
n'^fkij = n^nin"'njn^fkim + n^niT^vJ'fkim + rln"'njn'' fki^ + tIt^u^ fki^,. (4.37) 

In view of the boundary conditions (in the form (14.331) ). the two inner terms on the right- 
hand side of ()4.36|) and the two outer terms on the right-hand side of (|4.37j) vanish, and so 
the orthogonality is evident. 

Next we show that the boundary conditions are constraint-preserving. This is based on the 
following lemma. 

Lemma 25. Suppose that a and /?* vanish. Let Kij, f^ij he a solution to the homogeneous 
hyperbolic system l|3.36() . (|3.37() and suppose that the boundary conditions ()4.32|) are satisfied 
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on some face of dQ. Let Cj be defined by (|3.34|) . Then 

CjU^diC^ = (4.38) 

on the face. 

Proof. In fact we shall show that n^Cj = (so also n^Cj = 0) and r^n^diC^ = 0, which, by 
imphes (jl^ . First note that 

where we have used the first identity in H4.12() . Contracting with gives 

^iTi^i^k i I ^il^k^rri':^ . „k trim p 

where now we have used the equation (|3.37|) (with Pi = 0) for the first and last term and 
the second identity in (|4.12() for the middle term. From the boundary conditions we know 
that n^^n^n^ fkim = 0, and so the first term on the right-hand side vanishes. Similarly, we 
know that T^^n^Kim = on the boundary face, and so the second term vanishes as well 
(since the differential operator r^dk is purely tangential). Finally, n^S^'^fkim = n^{n^n"^ + 
pj^m _^ ^i^m-jj^.^ _ g^j^^ gQ ^Yie third term vanishes. We have established that n^Cj = 
holds on the face. 

To show that rJu^diC^ = on the face, we start with the identity 
Similarly 
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Therefore, 



For the last three terms, we again use (|ll(-i7|) to replace diKim with —fum and argue as before 
to see that these terms vanish. For the first term we notice that 5^^ dkdiKim = d^dkKim = 
Kim (from H3.36() and ()3.37|) with vanishing a and /?*). Since rP^n^Kim vanishes on the 
boundary, this term vanishes. Finally we recognize that the second term is the tangential 
Laplacian, r^^d^di applied to the quantity n^TP^Kira, which vanishes. This concludes the 
proof of (lOHl) . □ 

The next theorem asserts that the boundary conditions are constraint-preserving. 

Theorem 26. Let 0, be a polyhedral domain. Given initial data Kij{0), fkij{0) on 0, 
satisfying the constraints (|3.34j) and (|.S.40|) . define Kij and fkij for positive time by the 
evolution equations ()3.36|) . ()3.37|1 and the boundary conditions ()4.32() . Then the constraints 
(|3.34|) and (|3.4()j) are satisfied for all time. 

Proof. Exactly as for Theorem we find that Cj satisfies the wave equation (|3.4H) and 
both Cj and Cj vanish at the initial time. Define the usual energy 

E{t) = - [ {Cj& + d'CjdiC^) dx. 

Clearly E{0) = 0. From (|3.41j) and integration by parts 

E = [ CjU^diCUa. (4.39) 

Therefore, if a = and /3* = 0, we can invoke Lemma and conclude that E is constant 
in time. Hence E vanishes identically. Thus Cj is constant, and, since it vanishes at time 
0, it vanishes for all time. This establishes the theorem under the additional assumption 
that a and /?* vanish. 
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To extend to the case of general a and /3* we use Duhamel's principle. Let S{t) denote the 
solution operator associated to the homogeneous boundary value problem. That is, given 
functions Kjj(O), (j)kijiO) on Q, define S{t){Kij{0), (j)kij{0)) = (Kijit), cpkijit)) where 
is the solution to the homogeneous evolution equations 

^ij ^ 4'kiji 4' kij dk ^^ij i 

satisfying the boundary conditions and assuming the given initial values. Then Duhamel's 
principle represents the solution Kij, fi^-ij of the inhomogeneous initial-boundary value prob- 
lem (jOI)!) . (jOTl) . (lO^ as 

(Kijit), fkijit)) = S{t){Kij{0), A,,(0)) + r S{t - s){-d^dja{s),Lkijis)) ds. (4.40) 

Jo 

Now it is easy to check that the momentum constraint (|3.34|) is satisfied when Kij is replaced 
by didja{s) (for any smooth function a), and the constraint 1)3.40(1 is satisfied when f^ij is 
replaced by Lkij{s) defined by (|3.38|1 (for any smooth function /?*). Hence the integrand 
in (|4.4U|) satisfies the constraints by the result for the homogeneous case, as does the first 
term on the right-hand side, and thus the constraints are indeed satisfied by Kij, fkij. □ 

The analogue of Theorem El for the initial-boundary value problem follows from the pre- 
ceding theorem exactly as before. 

Theorem 27. Let 0, be a polyhedral domain. Suppose that initial data gij{0) and Kij{0) 
are given satisfying the Hamiltonian constraint and momentum constraint H3.34() . 

respectively, and that initial data fkij{0) is defined by (|3.39|) . Then the unique solution of 
the linearized EC initial-boundary value problem (|3.31() . ()3.36|) . (|3.37|) . together with the 
boundary conditions (|4.32|1 satisfies the linearized ADM system (|3.31|1 - (|3.34|1 in 

We close by noting a second set of boundary conditions which are maximal nonnegative 
and constraint-preserving. These are 

n'n^Kij = m'm^Kij = VP Kij = m'P Kij = n^n'm^ fkij = n'^n'V fkij = 0, (4.41) 
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or, equivalently, 

n'n^K.j = 0, T'W^'^K.j = 0, n'^nW^^fkij = 0. 

Now when we make an orthogonal expansion as in ()4.36|) . H4.37() . the outer terms on the 
right-hand side of the first equation and the inner terms on the right-hand side of the second 
equation vanish (it was the reverse before), so we again have the necessary orthogonality 
to demonstrate that the boundary conditions are maximal nonnegative. Similarly, to prove 
the analogue of Lemma I25( for these boundary conditions we show that the tangential 
component of Cj vanishes and the normal component of n^'diC^ vanishes (it was the reverse 
before). Otherwise the analysis is essentially the same as for the boundary conditions H4.32() . 

4.3.2 Extended EC System 

In this section we indicate an extended initial boundary value problem whose solution solves 
the linearized ADM system (|3.31j) - (|3.34|) in J7. This approach could present advantages from 
the numerical point of view since the momentum constraint is "built-in" , and so controlled 
for all time. The new system consists of H3.3ip . (|3.37p . and two new sets of equations 
corresponding to (|3.36jl 

Kij = -d^f^j + \ {diPj + djp^) - d^pAj - d^dja (4.42) 

and to a new three dimensional vector field pi defined by 

= d'Kii - diKl, (4.43) 

respectively. Observe that the additional terms that appear on the right-hand side of 1)4.42(1 
compared with (|3.36() are nothing but the negative components of the formal adjoint of the 
momentum constraint operator applied to pi. 

Let V be the vector space of quadruples of constant tensors {gij, Ky, fkij,Pk) symmetric 
with respect to the indices i and j. Thus diml/ = 33. The boundary operator An : V ^ V 
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in this case is given by 

Qij = 0, Kij = n^fkij - ^{niPj + UjPi) + n''pk6ij, fkij = rikKij, pi = -n^Ku + niK\. 

(4.44) 

The boundary operator An associated to the evolution equations (|4.42l) . (|3.37|) . and 

H4.43() has six positive, 21 zero, and six negative eigenvalues. Therefore, a nonnegative 
subspace is maximal nonnegative if and only if it has dimension 27. We claim that the 
following boundary conditions are maximal nonnegative for (|3.31j) . (|4.42j) . (|3.37)) . and (|4.43j) 

n'm^Kij = n'PKij = n^n^n^ fkij = n^{m'm^fkij +Pk) = n^{l'Vfkij +Pk) = n^mH^ fkij = 0. 

(4.45) 

These can be written as well: 

v}t^^K,,=Q, n'^n'n^ fkij = 0, n'=(r^V^-A,, + r'-p^) = 0, (4.46) 
and so do not depend on the choice of basis for the tangent space. 

Let us prove the claim that the subspace defined by (|4.45|) is maximal nonnegative. 
Obviously, dim = 27. Hence, it remains to be proven that N is also nonnegative. In view 
of (|4.44l) . the verification of nonnegativity of N reduces to showing that 

n^fkijK'^ - nVKi, + n'^pkKl > (4.47) 

whenever (|4.45() holds. In fact, we can prove that the left-hand side of 1)4. 47() vanishes 
pending 1)4. 45(1 holds. From the boundary conditions (in the form 1)4. 46() ') and the orthogonal 
expansions (|4.3(ij) and (|4.37j) of Kij and fkij, respectively, the first term on the right- 
hand side of (|4.47|) reduces to n^r'^W^'^ fkij Kim = —n^PkT^"^Kim- Then, combining the 
first and third terms of the left-hand side of ()4.47|) gives —n^pkT^-' Kij + n^pk5^^ Kij = 
n^PkU^n^ Kij. Finally, by using the orthogonal decomposition p> = n^pkU^ + T^-^Pk and the 
first part of the boundary conditions (|4.46|) the second term of the left-hand side of 1)4. 47() 
is —n^pkn^n^ Kij — pkrC'T^-^ Kij = —n^pkn^n^ Kij, which is precisely the negative sum of the 
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first and third terms of tlie left-hand side of (|4.47|1 . This concludes the proof of ()4.47|1 . 

Theorem 28. Let Q be a polyhedral domain. Suppose that the initial data gij{0) and Kij(0) 
are given satisfying the Hamiltonian (|3..S,3j) and momentum constraints (|3.34j) . respectively, 
fkij{0) is defined by H3.35() . and pi{0) = 0. Then the unique solution {gij,Kij, fj^ij^pi) 
of the initial boundary value problem (|3.31() . (|4.42l) . (|3.37|) . and ()4.43l) . together with the 
boundary conditions (|4.45j) . satisfies the properties pi = for all time, and (gij,Kij) solves 
the linearized ADM system ()3.31|) - (|3.34() in J7. 

Proof. Observe that the solution of the initial boundary value problem (|3.3ip . H3.36() . H3.37() . 
and (|4.32j) (boundary conditions), together with pi = for all time, is the unique solution 
of the initial boundary value problem H3.31|) . (|4.42j) . (|3.37p . and H4.43() . together with the 
boundary conditions H4.45() . The conclusion follows from Theorem 1271 □ 

We close by indicating a second set of maximal nonnegative boundary conditions (corre- 
sponding to dOB 'l for (jnSU, (11321), (ESZl), and dOH for which Theorem I2H1 holds as well. 
These are 

n^n^Kij = m^m^Kij = m^V K-ij = Pl^Kij = n^n^m^ fkij - rn^pk = n^n^V fkij - l^Pk = 0, 

(4.48) 

or, equivalently, 

nVK,, =0, r^'r^^K,,- = 0, nVr^Vfc^i " T^'Vfe = 0. (4.49) 
4.3.3 Inhomogeneous Boundary Conditions 

In this subsection we provide well-posed constraint-preserving inhomogeneous boundary 
conditions for (|3.31() . (|3.36|) . and ()3.37|) corresponding to the two sets of boundary condi- 
tions (|4.32|) . and (|4.41|) . respectively. The first set of inhomogeneous boundary conditions 
corresponds to (|4.32j) and can be written in the following form 



n'm^kij = n'Vkij = n^n'n^ fkij = n^m'm^ fkij = n'^l'P fkrj = n'^m'V fkij = 0, (4.50) 
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where Kij = Kij — Hij, fkij = fkij — Fkij^ with Kjj and F^ij given in i7 for all time and 
satisfying the constraints (|3.34|) and ()3.4U|) . respectively. The matter of choosing Kij and 
Fkij is deferred to the Appendix A. 

The analogue of Theorem 1261 for the inhomogeneous boundary conditions (|4.5Ujl is true. 

Theorem 29. Let Q be a polyhedral domain. Given initial data Kij{0), fkij{0) on $7 
satisfying the constraints (|3.,34|) and (|.S.4n|) . define Kij and fkij for positive time by the 
evolution equations ()3.36|) . (|3.37|) and the boundary conditions (|4.5U|) . Then the constraints 
H3.34() and (|3.4U|1 are satisfied for all time. 

Proof. Observe that Kij and fkij satisfy ()3.36|) and (|3.37|) with the forcing terms replaced 
by —didjU — d^'Fkij — kij and L^ij — dk^ij — Fj^ij, respectively. Exactly as in Theorem 1261 
it follows that Kij and fkij satisfy (|3.34|) and (|3.4()|) . respectively, for all time. Thus, Kij 
and fkij satisfy (|3.34|) and (|3.4fl|) . respectively, for all time. □ 

The analogue of Theorem for the case of the inhomogeneous boundary conditions 1)4.50(1 
follows from the preceding theorem by using the same arguments as in the proof of Theo- 
rem El 

Theorem 30. Let Q be a polyhedral domain. Suppose that initial data gij{0) and Kij{0) 
are given satisfying the Hamiltonian constraint (|3.33|) and momentum constraint (|3.34() . 
respectively, and that initial data fkiji^) is defined by (j3.39|l . Then the unique solution 
of the linearized EC initial-boundary value problem ((3.311) . ((3.36(1 . ((3.37(1 . together with 
the inhomogeneous boundary conditions ((4.50(1 satisfies the linearized ADM system ((3.31(1 - 
((!U?Ill in n. 

Note that there is a second set of inhomogeneous boundary conditions corresponding to 
((4.41(1 for which Theorem 1291 and Theorem 1301 remain valid. These are 

n'n^Kij = m'-m^kij = V-P Kij = m^H^Kij = n'^n'm^ fkij = n^n'-P fkij = 0, (4.51) 
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where again Kij = Kij — Kij , fkij = fkij — Fkij , with k^- and F^ij given and satisfying the 
constraints ()3.34|) and ()3.4U|) . respectively. 

Similar considerations can be made for the extended system introduced in the previous 
section. There are two sets of inhomogeneous boundary conditions for which the extended 
system produces solutions of the linearized ADM system (|3.31() - H3.34() on a polyhedral 
domain 17. These are 

n'm^Kij = n'Pkij = n^n'n^ fuj = n''{m'm^ fkij + Pk) = n^{l'V fkij + Pk) = n^m'V fkij = 

(4.52) 

and, respectively, 

n^n^Kij = m^m^Kij = m^V Kij = PPKij = n^n^m^ fkij - rri^pk = n^n^V fkij - l^Pk = 0, 

(4.53) 

where Kij and fkij are defined as before. 

The next theorem is an extension of Theorem |2H1 to the case of inhomogeneous boundary 
conditions. 

Theorem 31. Let Q be a polyhedral domain. Suppose that the initial data gij{0) and Kij{0) 
are given satisfying the Hamiltonian (|3.33() and momentum constraints ()3.34l) . respectively, 
fkij{0) is defined by (|3.35|) . and pi{0) = 0. Then the unique solution {gij, Kij, fkij, Pi) 
of the initial boundary value problem 1)3. 31() . ()4.42() . ()3.37|) . and ()4.43|) . together with the 
inhomogeneous boundary conditions (|4.52() ( or (|4.53l) ), satisfies the properties Pi = for all 
time, and {gij, Kij) solves the linearized ADM system (|3.31|) - (|3.34j) in Q,. 

Proof. Note that the solution of the initial boundary value problem 1)3. 31() . 1)3. 36() . 1)3. 37|) . 
and (|4.5()|) (or (|4.51|) . respectively), together with pi = for all time, is the unique so- 
lution of the initial boundary value problem (|3.3H1 . (|4.42jl . (|3.37|) . and (|4.43|) . together 
with the boundary conditions ()4.52l) (or 1)4. 53(1 . respectively). The conclusion follows from 
Theorem Oni □ 
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4.4 Alekseenko- Arnold Formulation 



4.4.1 Maximal Nonnegative Constraint-Preserving Boundary Conditions 

The equivalence of (|.S.50|) . (|.S.51|) and the hnearized ADM system has been studied in the 
second section of 3 for the case of pure initial value problem with the result that for given 
initial data satisfying the constraints, the unique solution of the linearized A A evolution 
equations satisfies the linearized ADM system, and so the linearized ADM system and 
the linearized A A system are equivalent (see Theorem 1.). Our interest is in the case 
when the spatial domain is bounded and appropriate boundary conditions are imposed. In 
this subsection, we provide maximal nonnegative boundary conditions for the linearized AA 
system which are constraint-preserving. This will then imply the analogue of the equivalence 
result proven in j^j (Theorem 1.) for the case of bounded domains. The ideas and techniques 
used in this section have many points in common to those used in Section 4.3 for the case 
of EC formulation. 

Assume that is a polyhedral domain. Consider an arbitrary face of d^l and let denote 
its exterior unit normal. Denote by m* and two additional vectors which together ri* form 
an orthonormal basis (see Figure Hm) . 

First we consider the following boundary conditions on the face: 

n^m^Kij = nH^Kij = n^ni^mP Xkij = ■nf'rPXkij = n^m^PXkij = 0. (4.54) 
These are equivalent to: 

n\^^Kij = 0, n^T'^T^'^Xkij = 0, (4.55) 

where tI := rriim^ + liV is the projection operator to n', and so do not depend on the 
choice of basis for the tangent space. We begin by showing that these boundary conditions 
are maximal nonnegative for the hyperbolic system (|3.5U() . 1)3. 51(1 . and so, according to the 
classical theory of jlH] and |^, the initial boundary value problem is well-posed. 
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Let V denote the vector space of pairs of constant tensors (Kjj, Xkij)- Here Kij is symmetric 
with respect to the indices i and j, and A/jj^ is a third-order constant tensor which is 
antisymmetric with respect to the first two indices and satisfies the cychc identity X^ij + 
Xjki + Xijk = 0. Thus dimy = 14. The boundary operator An associated to the evolution 
equations ()3.36|) . ()3.37|) is the symmetric hnear operator V ^ V given by 

kij = -n'^Afc(jj), Xkij = -n[kKiy. (4.56) 

A subspace of is cahed nonnegative for An if 

+ Xkij~X^'^ > (4.57) 

whenever {K,ij,Xkij) € N and {hij^Xkij) is defined by (|4.56j) . The subspace is maximal 
nonnegative if also no larger subspace has this property. Since An has five positive, four 
zero, and five negative eigenvalues, a nonnegative subspace is maximal nonnegative if and 
only if it has dimension nine. Our claim is that the subspace N defined by (|4.54|) is 
maximal nonnegative. The dimension is clearly nine. In view of (|4.56() . the verification of 
(|4.57() reduces to showing that n^XkijK^-' < whenever 1)4. 32() holds. In fact, vJ'XkijK^-' = 0, 
that is, n^Xkij and Kij are orthogonal (when (|4.54j) holds). To see this, we use orthogonal 
expansions of each based on the normal and tangential components: 

Kjj- = VL-mn^njKim + n'-riiTj^Kim + rln^njKim + rlrj^Kim, (4.58) 
n'^Xkij = vl-nivJ^njn^Xkini + n^riiT^n^Xkim + rlvJ^njU^Xkim + rlr^n'^Xkini- (4.59) 

In view of the boundary conditions (in the form (|4.55j) ). the two inner terms on the right- 
hand side of (|4.58j) and the last term of the right-hand side of (|4.59|) vanish. Also, the 
first two terms of the right-hand side of 1)4. 59() vanish due to the antisymmetry of Xkij with 
respect to the first two indices. Thus, the orthogonality is evident. 

Next we show that the boundary conditions are constraint-preserving. This is based on the 
following lemma. 
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Lemma 32. Suppose that a and /3* vanish. Let Kij, Xkij be a solution to the homogeneous 
hyperbolic system p.SUf) . (|3.51() and suppose that the boundary conditions ()4.54|) are satisfied 
on some face of dQ . Let Cj be defined by (|3.1fl|) . Then 

CjU^diC^ + n^Cjdi& = (4.60) 

on the face. 

Proof. In fact we shall show that n^Cj = (so also n^Cj = 0) and rJu^diC^ = 0, which, by 
imphes (jTBU]) . First note that 

Contracting with and using the first identity of ()4.12|) give 

where now we have used the equation H3.51() (with Pi = 0) for the last term. From the 
boundary conditions we know that n"^T^''dkKim = 0, and so the first term on the right- 
hand side vanishes (since the differential operator r'^^dk is purely tangential). Similarly, 
n^T^^Xkim = n'^T'^^'diKkrn = on the boundary face, and so the second term vanishes as 
well. We have established that n^Cj = holds on the face. 

To show that TjU^diC^ = on the face, we start with the identity 

TjH = r [n n + T )n = r n [o —Tj+rrn. 

Similarly 

80 



Therefore, 





The second term of the right-hand side vanishes since it involves the tangential Laplacian 
T^'^dkdi applied to the quantity TP"^n^Kijn, which vanishes. For the third and last terms, we 



resulting terms vanish. For the first term we notice that 6^^ dkdiKim = d^d^Kim = Him + 
d(iCm) — didmK\ (from (|3.5UI) and p.Slf) with vanishing a and /?*). So, r^'^n^S^^dkdiHim = 
rP'^n'kim + ^rP'^n'diCrn + ^rP'^n^dmCi - rP'^n'didmKl. From rP'^n'Kim = and n'd = 
on the face, this identity reduces to TP^n^5^^dkdiKim = ^TP™-n^diCm — TP™-n^didmn\- 

Finally, TP^nJ-n^n^dkdiHim = tP^{6^^ — T^'^)n"^dkdiKim = rP^n^dkd^ Kim — TP^T^^n^dkdiKim = 
TP^n^dkCm — TP^n^dkdmK\ — T^^T'-^n^^dkdiKim- Again, from n"^Cm = on the face and the 
tangentiality of rP^dk, the first term vanishes. Also, from the boundary conditions and the 
tangentiality of rP^r^'^dkdi, the last term vanishes. So, TP^v}v}n^dkdiK,im = —tP^ti^ dkdmi^\. 

From above, T^n'-diC^ = ^TP^n^diCm, which implies r^n^diC^ = on the face. This 
concludes the proof of (|4.60j) . □ 

The next theorem asserts that the boundary conditions are constraint-preserving. 

Theorem 33. Let Q be a polyhedral domain. Given 7ij(0), ACjj(O) on Q satisfying the 
constraints ^'^.'A'A\i and (jSISJ; define Xkij{0) by (|ll5l-i|) . Having Hij{0) and Xkij{0), define 
Kij and Xkij for positive time by the evolution equations (|3.5fl)) . (|.S.51|) . and the boundary 
conditions 1)4. 54|) . Then the constraints (|3.1U() are satisfied for all time. 

Proof. Exactly as in the proof of Theorem 1. of we find that Cj satisfies satisfies the 
elastic wave equation 



again use (|3.51|) to replace diKim with \/2Azim + diKim and argue as before to see that the 



Cj = d d(j_Cj) 



(4.61) 
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and both Cj and Cj vanish at the initial time. Define the usual energy 

E{t) = \ J^iCjC' + ^id'CjdiC^ + [di&f]} dx. 

Clearly E{0) = 0. From (|4.61|) and integration by parts 

E = [ {CydiC^ + n^Cjdi&)da. (4.62) 
Jan 

Therefore, if a = and /?* = 0, we can invoke Lemma |22 and conclude that E is constant 
in time. Hence E vanishes identically. Thus Cj is constant, and, since it vanishes at time 
0, it vanishes for all time. This establishes the theorem under the additional assumption 
that a and /?' vanish. 

To extend to the case of general a and /3* we use Duhamel's principle. Let S{t) denote the 
solution operator associated to the homogeneous boundary value problem. That is, given 
functions %(0), (ptijiO) on define S'(t)(%(0), 0fcij(O)) = (Oijit), cpkijit)) where (pkij 
is the solution to the homogeneous evolution equations 

-j=eij = d^(j)k{ij), -^i'kij = d[k0i]j, 

satisfying the boundary conditions and assuming the given initial values. Then Duhamel's 
principle represents the solution Kij, X^ij of the inhomogeneous initial boundary value prob- 
lem (j23ni), (ESII), (1^31) as 

1 

K(t),Afc»j(t)) = S(t)(K,j(0),Afeij(0))+ / Sit-s)i--^didja{s),r]kij{s))ds. (4.63) 







Now it is easy to check that the momentum constraint (|3.1U() is satisfied when Hij is replaced 
by —-^didja{s) (for any smooth function a), and, from (|3.50() . Cj{0) = \/2d^ {d^ r]i^(ij^ — 
^i^ki) ~ ^ ^'^y smooth shift vector Hence the integrand in 1)4. 63|) satisfies the 
constraints by the result for the homogeneous case, as does the first term on the right-hand 
side, and thus the constraints are indeed satisfied by Kj^, Xkij- □ 



82 



The analogue of Theorem 1. of |Hj for the initial boundary value problem follows from the 
preceding theorem exactly as there. 

Theorem 34. Let be a polyhedral domain. Suppose that initial data jij{0) and Kij{0) 
are given satisfying the Hamiltonian constraint (|3.9|) and momentum constraint (|3.1()|) . 
respectively, and that initial data Xkiji^) is defined by ()3.53l) . Then the unique solution 
of the linearized A A initial boundary value problem ()3.5U|) . ()3.51|) . and 7 defined by (|3.54() 
together with the boundary conditions (j4.54jl satisfies the linearized ADM system (j3.7|) - 

(|3Tnj) in n. 

We close by noting a second set of boundary conditions which are maximal nonnegative 
and constraint-preserving. These are 

n'n^Kij = m^m^Kij = rpKij = m^n'^n^Xkij = l^n^'n^Xkij = 0. (4.64) 

Now when we make an orthogonal expansion as in 1)4. 58(1 , 1)4. 59(1 , the first term on the right- 
hand side of the first equation and the first two terms on the right-hand side of the second 
equation vanish. The necessary orthogonality to demonstrate that the boundary conditions 
are maximal nonnegative follows from the boundary conditions and the antisymmetry of 
Xkij with respect to its first two indices. Similarly, to prove the analogue of Lemma IH^ for 
these boundary conditions, we show that the tangential component of Cj vanishes and the 
normal component of n^diC^ vanishes (it was the reverse before). The analysis is essentially 
the same as for the boundary conditions (|4.54j) . However, we need to do some more work 
here since the vanishing of both tangential component of Cj and normal component of 
v}diC^ is not enough for proving ()4.6U|) . We prove that djC^ vanishes on the face. From 
(|3.51|) . observe that X^^ = /\/2. So, from (|3.5()j) and the antisymmetry of Xkij in k and 
i, djC^ = V2djXf'^ = -V2djX^\ = -k\ = -5'^kij = - {n'n^ + m^m^ + PP)kij. From 
it follows that the divergence of vanishes on the face. Otherwise the approach follows 
exactly the ideas and techniques used for the boundary conditions (|4.54l) . 
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4.4.2 Extended A A System 



By using similar ideas and techniques as in Subsection 14.3.21 we indicate an extended initial 
boundary value problem corresponding to the AA formulation whose solution solves the 
linearized ADM system H3.7() - ()3.1U() in J7. The new system consists of (|3.7() for 7jj, p.Slfl . 
and two new sets of equations, one replacing (|3.5(J() 



and the other one corresponding to a new three dimensional vector field pi defined by 



Observe that the additional terms that appear on the right-hand side of H4.65() compared 
with (|3.5U|) are precisely the negative components of the formal adjoint of the momentum 
constraint operator 1)3. 1U|) applied to pi. 

Let V be the vector space of quadruples of constant tensors {'yij, Hij, Xkij,Pk), where 7jj, 
Hij are symmetric with respect to the indices i and j and X^ij belong to the 8-dimensional 
space T defined in Subsection 3.4.3. Thus dim = 23. The boundary operator An : V ^ V 
in this case is given by 

jij = 0, Kij = -n''Xk{ij)-^{niPj+njPi)+n'''pk6ij, Xkij = -n^j^K^^j, pi = -n'' Ku+niK,/ . 



The boundary operator An associated to the evolution equations (|3.7j) . (|4.(i5j) . (|3.51j) . and 
H4.66() has five positive, 13 zero, and five negative eigenvalues. Therefore, a nonnegative 
subspace is maximal nonnegative if and only if it has dimension 18. We claim that the 
following boundary conditions are maximal nonnegative for (|3.7|) , (|4.65|) , (|3.51|) , and (|4.66|1 



rim-^ Kij = rt'P Kij = {rri m-^ X^ij — Pk) = n'^^PPXkij — Pk) = n^rriVXkij = 0. (4.68) 




^didja, 



(4.65) 



Pi = d Kii - diKf. 



(4.66) 



(4.67) 
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These can be written as well: 



n'r^^Kij = 0, ^^(r^'r^^Afc,,- - r'^^) = 0, (4.69) 
and so do not depend on the choice of basis for the tangent space. 

Let us prove the claim that the subspace defined by (|4.68|) is maximal nonnegative. 
Obviously, dim = 18. Hence, it remains to be proven that N is also nonnegative. In view 
of ()4.67|1 , the verification of nonnegativity of N reduces to showing that 

Hijk'^ + Xhj~X'''' +Pif >0 (4.70) 

whenever 1)4. 68() holds. In fact, we can prove that the left-hand side of 1)4. 7U() vanishes 
pending that (|4.68|) holds: 

= —In'^K^^Xkij — In^p'Kij + 2n^piK\ 
= -2{n^n'n"'n^Kim + t^'t"^^ nim){rln"'njn^Xkim + T^T^n'^Xkim) - 2nY + 2n'pi4 
= —2n^T^^TpmKrpXkim — 2n^p' Kij + 2n^piK\ 
= -2n''pkT'^^Krp - 2n^p' Kij + 2n^piK\ 
= -2n''pk{S''"^ - ■rf'nF)Krp - 2nV'^ij + 2n^piK\ 
= 2n^ pkrf nP Hrp — 2n^p' Kij 
= 2n^pkinJ'vPKrp - 2-n}{n^pkn^ + T^^Pk)Kij 
= -2n^T^^ KijPk = 0. 

(4.71) 

This concludes the proof of 1)4. 7U() . 

Theorem 35. Let Q he a polyhedral domain. Suppose that the initial data jij{0) and Kij(O) 
are given satisfying the Hamiltonian H3.9|) and momentum constraints (|3.10() . respectively, 
XkijiO) is defined by ()3.53p . and pi{0) = 0. Then the unique solution {jijjKijjXkijjPi) 
of the initial boundary value problem (|3.7() . ()4.65|) . ()3.51|) . and 1)4. 66() . together with the 
boundary conditions (|4.68|) , satisfies the properties Pi = for all time, and {'jij , Kij ) solves 
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the linearized ADM system (|3.7() - p.lUf) in fi. 

Proof. Observe that the solution of the initial boundary value problem (|;i7|) , , (|H.51j) , 

and (|4.54j) (boundary conditions), together with pi = for all time, is the unique solution 
of the initial boundary value problem (|3.7() . (|4.65j) . H3.51|) . and H4.66() . together with the 
boundary conditions H4.68() . The conclusion follows from Theorem 1341 □ 

We close by indicating a second set of maximal nonnegative boundary conditions (corre- 
sponding to ((3211)) for (|T7|). (ESSl), (EISU), and for which Theorem 03 holds as well. 
These are 

n^npKij = m'm^Kij = rV Kij = m'^rin^Xkij + m^fc = l^rC'n^^kij + I'^Pk = 0. (4.72) 

4.5 Arnold Formulation 

Throughout this section, we will use the notations and results introduced in Subsection 3.4.4. 
In this part we unveil the technique that led us to the finding of maximal nonnegative 
constraint preserving boundary conditions for the three formulations, EC, AA, and A, 
respectively, analyzed in this thesis from this point of view. 

4.5.1 Maximal Nonnegative Constraint Preserving Boundary Conditions 

In this subsection we exhibit in detail our technique for finding maximal nonnegative con- 
straint preserving boundary conditions for Arnold's formulation. Everything that follows 
can be done for general polyhedral domains, but for the sake of simplicity let $7 be a par- 
allelepiped (a box) with faces parallel to the coordinate planes. Again, by m* and we 
denote two additional vectors which together n* form an orthonormal basis (see Figure 1^3)) . 

Denote hy v = Mu, and w = Mfi. Then, from (|3.61|) . we have that v and w are solutions 
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Figure 4.3: A parallelipipedic domain. 

of the following problem 

v = —- curlui, (4.73) 

w = - curliJ, (4.74) 

^;(0) = w{0) = 0. (4.75) 



Observe that (|4.73() . (|4.74() is a first order symmetric hyperbolic system. Its boundary 
matrix has the following eigenvalues (given with the corresponding eigenvectors): 

Ai = : (n, 0)^, (0, nf, 

A2 = -i : (m, 0^, (/_, -m)^, 

A3 = i : (1, mf i-m, l_f , 

where n is the exterior unit normal and m, I are chosen so that m _L n and I = n x m (see 
Figure 

If we consider maximal nonnegative boundary conditions for (|4.73|) - (|4.75|) . then the result- 
ing initial boundary value problem has the unique solution v = 0, w = 0. 
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From Proposition IIUI the general form of maximal nonnegative boundary conditions on the 
face having as exterior unit normal vector n = (1,0,0)"^ is given by 



(1 - I3)v2 + av3 + aw2 + (1 + f3)w3 = 0, (4.76) 
-6v2 + (1 + -f)v3 + (7 - l)w2 + 5w3 = 0, (4.77) 



with 



a (3 
7 5 



< 1. 



From 



V2 = 8x^12 + dy{-Vii - 7^33) + 9^2/23, 
^3 = dxVl?. + 5yJ^23 + dz{-Vll - 1^22), 
W2 = dxtJ-12 + dy{-Hu - m) + dzfl23, 

W3 = dxiJ-13 + dyii23 + dz{-mi - H22), 

the conditions (|4.76|1 - (|4.77|) become 

[(1 - /3)5^z^i2 + ad:riyi3 + aO^fiu + (1 + P)d^m]+ (4.78) 

dy[-{l - P)iyu - (1 - P)l^33 - afJ'33 + (1 + P)f^23] + 

dziil - (3)u23 - avii - OLV22 + afj.23 - (1 + (3)fJ'ii - (1 + (3)fJ-22] = 0, 

[-68^1^12 + (7 + l)d^i^i3 + (7 - 1)9x^12 + 6dccfii3]+ (4.79) 

dy[5uii + 51^33 + (7 + 1)J^23 + (1 - 7)/"ll + (1 - 7)^33 + '^>23] + 
dz[-6l^23 - (7 + 1)1^11 - (7 + 1)'^22 + (7 - 1)^23 - Sflll - 5fl22] = 

With no loss of generality (via Duhamel's principle), we can assume that the forcing terms 
in (HTTTll . inrTH]) vanish. 



From the system (|S?77|) - (|T75|) . it follows that 



5x1^13 = dyV2-i + dz{vii - 1^22) - 2/ii2, (4.80) 

d^Vi2 = dy{vii - 7^33) + a^z^23 + 2/ii3, (4.81) 

5x^13 = dy^i2z + djfiii - H22) + 2z>i2, (4.82) 

dxfJ-u = dyifiii - fi33) + dztJ-23 - 2z>i3- (4.83) 



By introducing (|I^ - (|X83|) into (|X78|) - (|I^ . we get 

dt[{l - P)m - am - avi3 + (1 + I3)ui2]+ (4.84) 



dy[{(5 - l)z^33 + ai/23 - a/"33 + (1 + /?)/^23] + 
5^[-ai/22 + (1 - ^)j^23 - (1 + /3)^22 + a/"23] = 0, 



dt[-5ixi3 - (7 + 1)mi2 - (7 - 1)^^13 + 5^12]+ (4.85) 

dy[8v3-i + (7 + 1)J^23 + (1 - 7)/^33 + <JAi23] + 
f^z[-(7 + 1)'^22 - 5V23 - 6fJ,22 - (1 - 7)^23] = 0, 

on the part of the boundary having (1,0, 0)^ as the unit exterior normal vector. 

Next, observe that we can get from the system p.77p - (|3.78p four relations that do not 
involve the transverse derivative dx 



dtvii + dy^ii3 - dziii2 = 0, (4.86) 

dt{i^22 + 1^33) - dyfiis + d,fii2 = 0, (4.87) 

dtm - dyvis + dzi^u = 0, (4.88) 

<9t(/i22 + /U33) + 5^1^13 - dziyu = 0. (4.89) 
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We will not use the relations H4.86|l and H4.88|) since they seem to introduce only complica- 
tions in what follows. 

By adding combinations of (|07|l and to both and we get 

dt[{l - /3)/xi3 - - on^n + (1 + P)i^i2 + w-z^22 + mi^ss + nii22 + nn33]+ (4.90) 

dy[{P - l)z^33 + - aiJ.33 + (1 + P)H23 " "1/^13 + ^1/13] + 

dz [-ai^22 + (1 - /3)l^23 - (1 + /?)At22 + a/"23 + "1^12 - = 



dy[6l^33 + (7 + 1)^^23 + (1 - 7)/^33 + 5/^23 " P^13 + ^2^13] + 
<92[-(7 + l)j^22 - SU23 - 6fJ.22 - (1 - 7)/^23 + ^2 " Q'J^12] = 

on the part of the boundary having (1,0, 0)^ as the unit exterior normal vector. 
It is easy to see now that (|4.84p - (|4.85p are satisfied if the following conditions hold: 



and 



dt[-SfJ,l3 - (7 + l)Ail2 - (7 - l)i^l3 + ^^12 +P1^22+P1^33 + 9/^22 + 9^33] + 



(4.91) 



mu22 + miy33 - avi3 + (1 + /?)z^i2 + n-^22 + n;U33 + (1 - ^3)^ll3 - a/Ui2 = 0, 



(4.92) 



(/? - l)z^33 + ah'23 + nui3 - a^33 + (1 + /3)/i23 - "2/U13 = 0, 



(4.93) 



az^22 + 1)J^23 + nui2 + (1 + /3)/U22 - am - "1/^12 = 0, 



(4.94) 



pU22 + P^33 + (1 - 7)Z^13 + S^12 + Q'Ai22 + Q'Ai33 " ^^13 - (7 + l)Ail2 = 0, 



(4.95) 



^^33 + (7 + l)i^23 + qi^l3 + (1 - 7)^33 + (5^23 - PAtl3 = 0, 



(4.96) 



(7 + 1)1/22 + 5z^23 + 91^12 + (^^^22 + (1 - 7)/^23 - PAil2 = 0. 



(4.97) 
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The matrix A corresponding to this hnear homogeneous system is: 



A 



m 


a 

P 


V 1 + 7 



m 

P-1 a 

(3-1 

p 

6 1 + 7 

6 



—a 


1 + /3 


n 


n 







1-/3 


—a 


n 








— a 


1 


+ /3 


— m 








n 


1 + /9 







-a 





— m 


1-7 


5 




q 







-5 


-1-7 


q 








1-7 




5 


-P 










5 





1 


-7 








Before we go any further, let us characterize the maximal nonnegative boundary conditions 
for (|H:77|) ~ (jT7H|) . 

The boundary matrix corresponding to this symmetric hyperbohc system has the following 
eigenvalues, given with the corresponding eigenvectors: 



A = 



n0n \ / \ / mQm+ l_ 
/ ' I nQn J \ 







m Q m + I Q I 







M 


f mQ n \ 




[ m Q n , 




[ -nQl 1 



A = 1 : 



^(/_0_/_ - mOm) 
m / 



m I 
\{m Q m - l_Q I) 





f mQ n 


M 


( n(dl_ \ 




[ l_Qn ) 




\ —m Q n 1 



-1 : 



m I 
\{lQl-rnQm) 



m / 
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From Proposition IIUI we know that if M G R^^^ and the norm of the matrix 



/ 1/^2 ^ 

1/^/2 

10 

\ 1 y 



M 



^ V2 ^ 

0^/2 

10 

1 y 





V 



is less or equal to 1, then the fohowing boundary conditions are maximal nonnegative: 



(m n)u + (_/_© njfi 
{n l)u — (m n)^ 
(m Du + _[ - m m)/2 
y u{m Qm-lQl)/2 + {rnQl)ix j 



\ ( 

+ M 



{l_ n)z/ + (jTT, n)/i 
(m n)i/ — (n 
1 - m m)/2 + (m _[)^ 



(4.98) 



For n = (1, 0, 0)^, m = (0, 1, 0)^, and / = (0, 0, 1)^, is equivalent to 



B 



( V22 ^ 



0, 



where 



B 



( -mi3 

-"T'23 
-"^33 
y 1 - m43 



mi3 2mi4 2mii 2(1 + m-12) 

772-23 2m24 2(1 + m2l) 2m22 

777-33 2(1 + 77134) 277231 277232 

-1+ 77243 277244 277741 277242 



77714 —"1-14 
-"1-24 
-1 + 77234 1 - 77234 



77244 



2?72i3 
277223 
277233 
2(1+77243) 



2(1 - 77212) 277211 

-277222 2(-l + 77221^ 
-277232 2772-31 
— 277242 277241 



\ 



-77244 

(4.99) 



Now, we want a match between the systems ()4.92() - ()4.97() and 1)4. 98() . First of all, observe 
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that the first two columns of B coincide except the sign and, since we want each row of A 
to be a hnear combination of rows of B, this induces 

a = 0; /3 = 1; 7 = -1; 5 = 0; m = 0; p = 0. 

Also, we have to worry about the rank of the matrix A, which is a 6 x 10-matrix. Obviously, 
the rank of A cannot exceed four. Fortunately, we know that the rank of A equals the rank 
of AA^ which is a 6 X 6-matrix. Using this fact, we can see that A has rank four for the 
above values of the other parameters if and only ifn^ + q^ = 4 (otherwise A has rank greater 
or equal to five). Prom here, we get that the rows of B must be linear combinations of the 
rows of A too and it turns out that M must be 



M = 



/ 





1 








\ 




1 
























-1 




V 








1 





/ 



So, the maximal nonnegative boundary conditions on the face are: 

1^12 = 0, (4.100) 

uis = 0, (4.101) 

-M22 + M33 = 0, (4.102) 

/X23 = 0. (4.103) 

Same kind of conditions could be obtained for each face of fi. In fact, we can write the 
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boundary conditions for all faces in the following form: 

= 0, (4.104) 

= 0, (4.105) 

^ [ 0, (4.106) 

iQl— mQml \ jj, 

\ f u 



nQm 






























mQ I I \ fi 



0, (4.107) 



or 



e X, (4.108) 



where 



X = Span{ 



= / \dn 



nQn \ I mQl_ \ I m Q m \ I m Q m — l^Q l_ \ I 
/'\ /'\ /'\ I \ nQn 



\ \ 

nQm I \ nQl^ I \ rnQm + l_Ql_ 

Observe that X is a maximal nonnegative (in fact null) subspace for the boundary operator 



An 



Kn 

-K„. 



where by definition 

KnU = -[Skw(n)u — uSkw(n)], 



= T = 



for all u e Mf^^. 
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Also, it is easy to see that the boundary conditions H4.1U4|) - (|4.1U7|1 can be written as: 



nW^'^Uij = 0, (2t^V^' - T''^T'^)fiij = 0, (4.109) 

where r*-' := m^m^ + is the projection operator orthogonal to n', and so, there is no 
dependence on the particular choice of m* and /*. 

Define 

E{t) = / lup + \w\'^dx. 
Jn 

Then, 

E{t) = / w{n X v)da = — v{n x w)da, 
Jdn ~ ~ ~ J an 

where n is the exterior unit normal vector. 

By using (|4.104|) - (|4.1fl7j) . it is easy to prove that E{t) = and this implies again that 
Mu = and My, = for all time. 

By the stucture of the main system of differential equations ()3.77|) - ()3.78|) . there is a second 
set of maximal nonnegative constraint preserving boundary conditions 

nV^ Vij = 0, (2t^V^' - r'^'r^^Vij- = 0. (4.110) 

Moreover, by our above analysis there are no other sets of maximal nonnegative constraint 
preserving boundary conditions corresponding to the system H3.77() - ()3.78|) . 

4.5.2 Equivalent Unconstrained Initial Boundary Value Problem 

Consider the constrained initial boundary valu problem i|2.2()j) - (|2.22)) together with maximal 
nonnegative boundary conditions 

uian e Y, (4.111) 
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where y is a maximal nonnegative boundary space. Also, consider the associated extended 
problem 1)2. 23() - (|2. 24(1 with boundary conditions 

dY = Y x{BnY)^, (4.112) 

\dn 

where Bn{x) = —Ti^^^B^rij is the boundary matrix corresponding to the first order differ- 
ential operator B at x ^ dO,. 

From Corollary El ^ is nonnegative for the boundary matrix A„ associated to ()2.23() if and 
only if Y is nonnegative for the boundary matrix An associated to ()2.2U() . 

Furthermore, assume Y is maximal nonnegative for A^. Let t; be a solution of 1)2. 20[) - (|2. 22(1 
and (|4.111l) . and {u,z)'^ be a solution of (IT^ - dT^ and (14.1121) . 

A formal calculation shows that 

^^(Ik - v\\h(n) + \\4h(n)) = -\ J^^iu - vfAniu -v)da < 0. 




From here, it follows that z = and u = v. In other words, the initial boundary value 
problems ^T^-^TT^ with (14.1111) and (tT^ - (IT^ with MlVll are equivalent. 

Returning to our problem, observe that the boundary matrix of the extended system ()3.90p 
reads 



where 




with MnU = u n — (tiu) n, Vn S M^, 
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Straighforward computations show that 



BnX = Span{ ( ^ ) ' ( ) ' ( ^ ) ' ( " I ^' 



and so, 



{Br,X)^ = Span{ ( - j ; |^ ° | }. 

Therefore, X = X x (BnX)-^ has dimension ten and it is nonnegative for A^. Since An 
has exactly ten nonnegative and eight strictly negative eigenvalues, it follows that X is 
also maximal nonnegative. Thus, the following boundary conditions are suitable for the 
extended system (|8.9n|) 



n m 

/ V 

I \ fi 


I (d I — m(d m j ^/i 




= 0, (4.113) 
0, (4.114) 
= 0, (4.115) 

(4.116) 



m Q I I \ fi 

p = 0, (4.117) 
qn = 0, (4.118) 



or, 



nV-'Vij = 0, (2t^V^' - T^W'^)nij =0, pi = 0, n'qi = 0. (4.119) 

Note the existence of a second set of boundary conditions for the extended system (|3.9Up 
related to (|4.11()|) : 

nV^ Vij = 0, (2r*V^'' - T^^T'^)vij = 0, pi = 0, n'qi = 0. (4.120) 
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Appendix A 



General Solution for the Linearized 
Momentum Constraints 



In order to apply the considerations of Subsection 4.2.3 to the EC problem, it would be 
useful to solve the linearized momentum constraint equations ()3.34|) for symmetric matrix 
fields defined on the entire space. 



Define the symbol ejjfc by 



Note that 



+1 if (ijk) is either (123), (231), or (321), 

-1 if (ijk) is either (321), (213), or (132), (A.l) 

otherwise. 



From Cj{K) := d^Kji — djK\ = d^Kji — K'lSji) = 0, there exists a matrix field aij so that 
Kjk — K\6jk = ^k^diaji. Therefore, 

= e/9iCTj7 - ^eP'^diapi6jk. (A.3) 
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Since Kij is symmetric, 

eC^diaji = e.diaki. (A.4) 



It is not difficult to see that HA.4|1 is equivalent to 

d^mi - ^^al = 0. (A.5) 

The proof of this fact follows easily by multiplying ()A.4|) with e'^-'*, and using the identity 
(US). 

From ()A.5|) and d^au — dial = d\aii — (tISu), there exists a matrix field rjij such that 

(Tij - cfl^ij = ei^^dkmj, (A.6) 

and so, 

cTij = el'^dkVij - ^e"''^dkVisSij. (A.7) 

By substituting ()A.7|) into (|A.3|) we get the general solution of the linearized momentum 
constraint H3.34|) . 
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Appendix B 

The 4— D Covariant Formulation 



Assume that the metric g^i^ is a shght perturbation of the Minkowski metric ijni,: 

Qixv = ViJ.u + h^^, \h^^\ < 1. 

To first order in /i, the Christoffel symbols are given by 



and so, to first order in h, the Ricci tensor is 
where □ is the d'Alembertian. 

Therefore, the Einstein equations in first approximation read 

- - Km + ^A,M- = -le^rGV, (B.3) 

where 

S^v = T^i, — —rini/T^ 
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is the source term. 



Here T^jy is taken to lowest order in /i^j, and satisfies the ordinary conservation conditions 

r,^M = 0. 



We can not expect a field equation such as HB.3|) to provide unique solutions because given 
any solution, we can always generate other solutions by performing coordinate transforma- 
tions. The most general coordinate transformation that leaves the field weak is of the form 

x^' — .x'^ = x'^ + e^(x), (B.4) 
where e(t is at most of the same order of magnitude as /i^jy . 
From 

gt^U = l^t^'^ - /iM^'^ (B.6) 



it follows that 



where = 'q^^e" . 



hflU hflU ^fJ-,!^ ^I'l/i' (-'^•^) 



By direct computations, we can prove that HB.8|) is also a solution of ()B.3|) . To remove this 
non-uniqueness of solution, we have to work in some particular gauge. It turns out that 
one of the most convenient choice is the harmonic gauge, for which 



= g^'^r^^u = 0. (B.9) 



Let us prove that there is always a coordinate system in which ()B.9|) holds 
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Recall that the Christoffel symbols transforms as 



_ dx'^ dx^ dx" dxP dx" d^x'^ 

9x^9x^5^ ~ dx'^ dx't' dxPdx'' ' ^ 



Contracting this with g''^'^ , we get 



dx'^ d'^x'^ 



Thus, if does not vanish, we can always find a new coordinate system x"^ by solving the 
second-order partial differential equations 

dx'^„„ _ d^x'^ 



dxP^ ^ dxPdx'' ^' (B.12) 
which gives T'^ = in the x'-coordinate system. 

So, without loss of generality, assume that the chosen gauge is harmonic. Then, since 
= 0, it follows that the coordinates are harmonic functions 

□x^ = {g^''x''.).k = 9'''' fll\ - = -TP = 0. (B.13) 

This explains the term harmonic gauge. 
For a harmonic gauge, the equation (|B.3|) is 



□V = -IQttGS^^. (B.14) 



One solution of this equation is the retarded potential 

h,,{x,t) =4G [ M^^liZ-^^dx'. (B.15) 
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To this solution, we can add any solution of the corresponding homogeneous equation 



which also satisfies 

d Id 

K = -^ir:zK. (B.17) 



which comes from = to first order. 

Observe that Hij.l5|) automaticahy satisfies the harmonic condition (|B.17() (at least for a 
compactly supported source S^i^"), since 



that follows from the conservation law 

d 



T!f = 



and 

Interpretation: The solution ()B.15|) is interpreted as the gravitational radiation produced 
by the source S^^u, whereas any additional term (solution of HB.16|) . and ()B.17|) ) represents 
the gravitational radiation coming from infinity. The term t — \x — x'\ in (|B.15() shows that 
the gravitational effects propagate with the speed of light. 

The drawback of this method is the fact that the gauge condition T'^ = can be imposed 
only locally 
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